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Abstract 



In this thesis we investigate cosmological models more general than the 
isotropic and homogeneous Friedmann-Lemaitre models. We focus on cos- 
mologies with one spatial degree of freedom, whose matter content consists of 
a perfect fluid and the cosmological constant. We formulate the Einstein field 
equations as a system of quasilinear first order partial differential equations, 
using scale-invariant variables. 

The primary goal is to study the dynamics in the two asymptotic regimes, 
i.e. near the initial singularity and at late times. We highhght the role of spa- 
tially homogeneous dynamics as the background dynamics, and analyze the 
inhomogeneous aspect of the dynamics. We perform a variety of numerical 
simulations to support our analysis and to explore new phenomena. 
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Chapter 1 
Introduction 



1.1 Background in cosmology 

The goal of relativistic cosmology is to describe the structure of the observ- 
able universe on a sufficiently large scale, and to explain its origins, using 
Einstein's General Theory of Relativity. 

The most fundamental cosmological observation is that light from dis- 
tant galaxies is redshifted (this was discovered by Hubble in the 1920s). The 
standard interpretation of the redshifts is that the universe is expanding. 
The second fundamental observation is that the universe is filled with radi- 
ation (as discovered by Penzias & Wilson in 1965), which has a blackbody 
spectrum to a very high accuracy, corresponding to a present temperature 
of about 3 K. The standard interpretation is that this radiation, called the 
cosmic microwave background (CMB), was emitted when the universe was 
denser and hotter, and has travelled freely to us ever since. A third obser- 
vation, that provides information about the physics of the universe at early 
times, concerns the overall abundance of helium in the universe (see Pee- 
bles 1993 inOl pages 143-4] and Hoyle & Tayler 1964 The constant 
abundance of helium in stars and intergalactic matter (^ 25% by mass) over 
time suggests that helium was formed during a very early epoch when the 
universe was very hot. These three observations form the physical basis for 
the so-called hot big-bang theory of the universe, which can be summarized 
in the following history of the universe. 

The history begins at the so-called Planck time (see, for example. Coles 
& Lucchin 1995 page 111]), the time after which classical general rela- 
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tivity may be assumed to be valid. The corresponding Planck temperature 
is approximately 10^^ K. At such a high temperature, photons have enough 
energy to produce matter and anti-matter in collisions. When the tempera- 
ture fell to the order of 10^ K, collisions of photons could no longer produce 
even the lightest known particles and antiparticles (electrons and positrons), 
and the radiation-dominated epoch began (i.e. most of the energy was in 
the form of radiation). Slight asymmetries between matter and anti-matter 
resulted in all anti-matter being annihilated by matter, and the remaining 
matter consisted of protons, neutrons, electrons, photons and neutrinos. It 
was also cool enough at this temperature for most neutrons to be rapidly tied 
to form helium nuclei in a process called nucleosynthesis (see, for example, 
Weinberg 1977 ^3 page HO]). Radiation was "coupled" to matter in the 
sense that photons were frequently scattered by electrons and nuclei through 
Compton scattering, and matter was kept ionized. As the universe expanded 
and cooled, radiation lost energy faster than matter did, with the result that 
eventually the universe became matter-dominated. When the temperature 
dropped to about 3000 K, photons no longer had enough energy to con- 
stantly ionize matter. In other words, the free electrons became bound to 
nuclei, with the result that the photons could travel freely. One says that 
the photons have decoupled from the matter. These photons form the CMB 
radiation, which is detected at the present time. The CMB thus provides 
a snapshot of the universe at a temperature of about 3000 K, but it is also 
affected by the geometry of spacetime as it travels towards us. 

Within this framework, the formation of galaxies occurs due to the growth 
of density fluctuations under the influence of gravity. The key point is that 
random fluctuations in the matter density do not grow rapidly enough to 
form galaxies (see, for example, Peebles 1993 jHOl page 619]). Hence one 
has to postulate the existence of a suitable spectrum of density fluctuations 
at very early times. The most popular, but not universally accepted, ex- 
planation of the origin of density fluctuations is the theory of inflation in 
the early universe, first proposed by Guth 1981 During the inflationary 
epoch, which is postulated to have occurred shortly after the Planck time, the 
universe underwent exponential expansion, and quantum fluctuations were 
hugely magnified to create a spectrum of density fluctuations, which were 
sufficiently large to form galaxies, but not so large as to produce too many 
black holes (see also Peacock 1999 [59, Chapter 11]). 

Over the past 20 years, as cosmological observations have become more 
extensive and more accurate, an increasingly complicated picture of the uni- 
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verse has emerged. Redshift surveys (e.g. the Sloan Digital Sky Survey jHS], 
the Las Campanas Redshift Survey jlHl) show that galaxies are not dis- 
tributed uniformly, but form clusters which in turn are grouped into huge 
structures on a scale of 100 Mpc^, e.g. the Great Attractor (probably a super- 
cluster, with galaxy cluster Abell 3627 near its centre), which is attracting 
several other clusters, including the Local Group. These observations are 
being used to determine the power spectrum of the distribution of galaxies 
(see Strauss 2004 |68j, Tegmark et al. 2004 [69J) and to compare it with 
the CMB power spectrum (see Hu & White 2004 [40 ). Until recently, it 
was widely believed that the rate of expansion of the universe is slowing 
down. This belief was challenged when the recent observations of type la 
supernovae suggested an accelerating expansion (see, for example, Riess & 
Turner 2004 [63^ and Turner & Riess 2002 |71j). This observation, when com- 
bined with the observations of the power spectrum of the CMB temperature, 
restricts the total density parameter^ Vt to have a value close to 1 (see, for 
example, Spergel et al. 2003 jSHI)- However visible matter contributes a value 
of approximately 0.005 to f2, and non-luminous baryonic matter contributes 
a value of approximately 0.035 [5H]. This dramatic conclusion suggests that 
the matter-energy content of the universe is dominated by yet unidentified 
non-baryonic dark matter and/or energy, indicating that our understanding 
of the physical universe is still quite incomplete. 

The simplest mathematical models within the framework of classical gen- 
eral relativity that describe an expanding universe are the Friedmann-Lemaitre 
(FL) models^, which are spatially homogeneous and isotropic (see Ellis 1989 
for a discussion of the historical development). An important question is this: 
Do cosmological observations support the assumption of isotropy and homo- 
geneity on a sufficiently large scale? In view of the innate limitations in 
making cosmological observations (we can observe only from one point in 
space and one instant in time on a cosmological scale), and the difficulties in 
interpreting the observations, the answer is not clear-cut. It is thus desirable 
to investigate the dynamics of cosmological models more general than the 
FL models, and to determine what restrictions observations impose on them. 
This general area of investigation forms the basis of this thesis. 



^1 parsec is about 3 lightyears; the size of a typical galaxy is 1 Mpc and of a typical 
cluster is 10 Mpc. 

^The density parameter is defined in Section H"^ 
■^These models are introduced in Section [TT^ 
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1.2 Cosmo logical Hierarchy 

In this section, we define cosmological models and classify them according to 
their symmetries. 

Cosmological models 

A cosmological model (A^,g,u) is determined by a Lorentzian metric g de- 
fined on a manifold A4, and a family of fundamental observers, whose con- 
gruence of world-lines is represented by the unit timclike vector field u, which 
we often identify with the matter 4- velocity. The dynamics of the model is 
governed by the Einstein field equations (EFEs) 

Gab — Tab, ^ (1-1) 

with suitable matter content. 

The matter content is assumed to be a perfect fluid with energy density 
> 0, pressure p and (unit timelike) fluid 4- vector u: 

Tab = I^UaUb+p{gab + UaUb). (1.2) 

We assume that the equation of state of the perfect fluid is of the form 
p = (7 — 1)|U, where 1 < 7 < 2 is constant. The cases 7 = 1 (dust) and 7 = | 
(radiation) are of primary physical interest. The cosmological constant A in 
the EFEs can be treated as a perfect fluid with /x = A and p = —A. A more 
formal introduction will be presented in the next chapter. 

It is helpful to classify cosmological solutions of the EFEs using the di- 
mension of orbits of the symmetry group admitted by the metric. This 
classification scheme forms a hierarchy of cosmological models of increasing 
complexity, starting with the FL models: 

1) Friedmann-Lemaitre (FL) cosmologies {Gq) 

2) Spatially homogeneous (SH) cosmologies (G3) 

3) G2 cosmologies 



^We use geometrizcd units throughout, i.e. c = 1 and SwG = 1, where c is the speed 
of hght in vacuum and G is the graviational constant. 
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4) generic cosmologies (Go) 

An FL cosmology admits a local 6-parameter group of isometrics acting 
on spacelike 3-surfaces. i.e. it has 6 independent Killing vector fields (KVFs) 
whose orbits form spacelike 3-surfaces. An SH cosmology is either a Bianchi 
cosmology, which admits a local 3-parameter group of isometrics acting 
simply transitively on spacelike 3-surfaces, or a Kantowski-Sachs cosmology, 
which admits a group G4 of isometrics acting on spacelike 3-surfaces but 
does not admit a simply transitive G3 subgroup. A G2 cosmology admits a 
local two-parameter Abelian group of isometrics with spacelike orbits, per- 
mitting one degree of freedom as regards inhomogeneity. From the viewpoint 
of analysis, the equations for G3 cosmologies are ordinary differential equa- 
tions (ODEs), those for G2 cosmologies are 1-1-1 partial differential equations 
(PDEs), and those for Go cosmologies are 1-1-3 PDEs. 



Friedmann-Lemaitre cosmologies 

As mentioned in the introduction, FL models are spatially homogeneous 
and isotropic. In physical terms, "spatially homogeneous" means "the same 
at every point at any given instant of time", while "isotropic" means "the 
same in every spatial direction" . Mathematically, the requirements of spatial 
homogeneity and isotropy imply that the line element admits the maximal 
group Gq of isometrics acting on spacelike 3-surfaces. With an appropriate 
choice of time coordinate, the line element can be written as ^ 

ds^ = -dt^ + l{tfga/3dx''dxf^, (1.3) 

where ga/sdx'^dx^ is the metric of a (spacelike) 3-surface with constant curva- 
ture (see, for example, Weinberg 1972 [H21 Section 13.5]). The function l(t) 
is called the length scale of the model. 

The two most important quantities describing an FL model are the Hub- 
ble scalar H, where H = I /I, and the total^ density parameter, defined by 

n = -^. (1.4) 



■^Latin indices run from to 3; Greek indices from 1 to 3; capital Latin indices run 
from 2 to 3. 

^In general, /x is a sum of terms, for example matter, radiation and the energy density 
due to a cosmological constant. 
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FL models are further separated into three types according to the curvature 
of the spacehke 3-surfaces t = const., which is completely determined by Q: 

Open FL: negative curvature, f2 < 1 

Flat FL: zero curvature, Q = 1 

Closed FL: positive curvature, f2 > 1 

We want to be able to determine whether a given cosmological model 
is exactly FL or close to FL. For this purpose it is necessary to introduce 
physical quantities called the kinematical quantities (see, for example, Wain- 
wright & Elhs 1997 (WE) [771 page 18], and Ellis 1971 [HI page 111]). For a 
unit timelike vector u, one can decompose the covariant derivative Ua-^b of the 
corresponding 1-form into its irreducible parts - the kinematical quantities 
o'ab, ^ab, © and iia- The scalar 6 is the rate of expansion scalar and iia is the 
acceleration vector. The tensor aab is the rate of shear tensor and ujab is the 
rate of vorticity tensor. If the vorticity is zero, we say that u is irrotational. 
Because we are working in a cosmological context we shall usually replace 
O by the Hubble scalar H defined earlier, which is related to B according 
to if = |0. The following theorem characterizes the FL models using the 
kinematical quantities (see, for example, Krasihski 1997 IIH page 11]): 

Theorem 1 

A solution of the EFEs with perfect fluid matter content is an FL 
model if and only if 

C^o6 = = UJab = Ua- (1-5) 

The FL models can also be characterized using the Weyl conformal cur- 
vature tensor Cabcd, which is interpreted as describing the free gravitational 
field. The Weyl tensor can be decomposed into its electric part Eab and 
magnetic part Hab relative to u. The following theorem characterizes the FL 
models using Eab and Hab (see, for example, Wainwright 1996 [7^1 Section 
4.1]): 

Theorem 2 

A solution of the EFEs with perfect fluid matter content and equation 
of state p = p{fi) is an FL model if and only if 



Eab = = Hab- 



(1.6) 
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Spatially homogeneous cosmologies 

Spatially homogeneous cosmologies with an isometry group G3 acting on 
spacelike 3-surfaces, the so-called Bianchi cosmologies, are classified accord- 
ing to the structure of the G3 group, as given in Table ITTTl ^ 

Table 1.1: Classification of Bianchi cosmologies. 



Group type Eigenvalues 
Class A Class B of Uajj 



IX 




+ 


+ 


+ 


VIII 




+ 


+ 




VIIo 


VII, 


+ 


+ 





VIo 


VI, 


+ 







II 


IV 


+ 








I 


V 












Bianchi cosmologies comprise two classes, 

i) the non-tilted models in which the 4- velocity vector u is always orthog- 
onal to the spacelike G3 orbits, were first studied in detail by Ellis & 
MacCallum 1969 23 . 

ii) the tilted models in which u is never orthogonal to the spacelike G3 
orbits, were first studied by King & Ellis 1973 |I3] . 

The dynamics of non-tilted Bianchi cosmologies have been analyzed in 
detail (see WE, Hewitt et al. 2003 |36, and Horwood et al. 2003 [38|). Much 
less is known about the dynamics of tilted Bianchi cosmologies and work is 
currently in progress. * Indeed in this thesis we provide a formulation of 
the EFEs for analyzing the asymptotic and intermediate dynamics of the 
so-called G2-compatible Bianchi cosmologies. 

There exist G3 subgroups within the Gq isometry group of FL cosmologies. 
Thus FL cosmologies appear as special cases of SH cosmologies - open FL as 

^The symmetric matrix in Table [TTTI is defined in equation (|2.3|l . 
®See Hervik 2004 j321 and Hervik & Coley 2004 
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Bianchi V and VII/i, flat FL as Bianchi I and VIIq, and closed FL as Bianchi 
IX (see WE, page 37). Notice that open and flat FL cosmologies appear as 
more than one type of Bianchi cosmology. These are instances of multiple 
representation of special models in a more general class of models. 

Spatially inhomogeneous cosmologies 

The simplest spatially inhomogeneous cosmologies are the G2 cosmologies, 
which admit one spatial degree of freedom. These models, which form the 
main focus of this thesis, will be classifled in Chapter IHl We do not consider 
generic spatially inhomogeneous cosmologies, i.e. the so-called Go cosmolo- 
gies, in this thesis. We refer to Uggla et al. 2003 |72] and Lim et al. 2004 UHl 
for a general framework for analyzing these models. 

1.3 Goals 

We now discuss three physically important goals of the study of cosmological 
models more general than FL. 

PI. Describe the asymptotic dynamics at early times. 

Assuming that the universe is expanding from an initial singularity, 
the phrase "dynamics at early times" refers to the dynamics in the 
asymptotic regime associated with the initial singularity, characterized 
by / — >^ 0, where / is the overall length scale. Classical general relativity 
breaks down at the singularity, and should be replaced by a theory of 
quantum gravity there. The time after which classical general relativity 
is valid is called the Planck time. Thus in physical terms the phrase 
"dynamics at early times" refers to the evolution of the universe in an 
epoch starting at the Planck time. It is important to classify all possi- 
ble asymptotic states at early times that are permitted by the EFEs, in 
order to explain how the universe may have evolved. 

Two important conjectures about the asymptotic dynamics at early 
times have been made by Belinskii, Khalatnikov & Lifshitz, which we 
shall refer to as BKL I and BKL II. 
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BKL conjecture I 

For a typical cosmological model the matter content is not 
dynamically significant near the initial singularity (Lifshitz 
& Khalatnikov 1963 page 200], BKL 1970 [HI pages 
532, 538]). 

The implication of this conjecture is that one can use the vacuum EFEs 
to approximate the evolution of a typical cosmological model near the 
initial singularity. An example of atypical cosmologies is the set of mod- 
els with an isotropic singularity (Goode & Wainwright 1985 [22], Coley 
et al. 2004 ^j). Another example is the set of models in which the 
matter content is a so-called stiff fluid {p = /i) (Andersson & Rendall 
2001 HI). 

BKL conjecture II 

Near the initial singularity, the spatial derivatives in the 
EFEs for a typical cosmological model are not dynamically 
significant (BKL 1982 page 656]). 

The implication of this conjecture is that near the initial singularity the 
EFEs effectively reduce to ODEs, namely, the ODEs that describe the 
evolution of SH cosmologies. Hence it is plausible that SH models will 
play an important role in analyzing spatially inhomogeneous models. 

P2. Describe the asymptotic dynamics at late times. 

The phrase "dynamics at late times" refers to the dynamics in the 
asymptotic regime characterized by / — oo. For example, the dynam- 
ics at late times in a radiation model would describe the later stage 
of the radiation-dominated epoch before the universe became matter- 
dominated. Likewise, the dynamics at late times in an inflationary model 
would describe the final stage of the inflationary epoch. 

P3. Describe the dynamics of close-to-FL models. 

As mentioned earlier, it is usually assumed that deviations from exact 
isotropy and spatial inhomogeneity can be described using linear per- 
turbations of an FL model. ^ Linear perturbations have been used in 

^For a review of cosmological perturbations, see Mukhanov et al. 1992 Other 
references that use a gauge-invariant approach are Bardeen 1980 W , Bruni et al. 1992 Jl] 
and Goode 1989 
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studying, for example, the evolution of the primordial density pertur- 
bations during the radiation-dominated epoch, and the development of 
large-scale structure in the distribution of galaxies. A key question is 
this: does linear perturbation theory correctly describe the dynamics of 
close-to-FL models? 

1.4 Overview 

This thesis is divided into two parts. The first part consists of Chapters |2l- 
13 which provide a unified formulation of the EFEs, the background and 
analytical results for the dynamics of G2 and SH cosmologies. The second 
part consists of Chapters ISHHl in which we present a numerical exploration 
of the dynamics of G2 cosmologies. 

We first establish a framework for the study of inhomogeneous cosmolo- 
gies in Chapter |21 We use the orthonormal frame approach to formulate the 
EFEs, and introduce the appropriate normalizations for G2 and SH models 
respectively, which lead to scale-invariant variables that are bounded at the 
initial singularity. 

In Chapter 121 we present a hierarchy of subclasses within the G2 cosmolo- 
gies, and develop the links between G2 and SH cosmologies. 

A new, unified formulation for G2-compatible SH cosmologies is given in 
Chapter m We discuss the tilted SH cosmologies in detail, emphasizing the 
past asymptotic dynamics and giving conditions for the existence of Mixmas- 
ter dynamics. 

Chapters |SHH1 deal with analytical aspects of G2 cosmologies. In Chap- 
terlSlwe complete the formulation for G2 cosmologies by choosing appropriate 
temporal and spatial gauges. We discuss the notions of the silent boundary 
and asymptotic silence, which play a key role in linking G2 dynamics with 
SH dynamics. Chapter IHl presents two G2 solutions - one of which is new 
- that develop large spatial gradients as they approach the initial singular- 
ity. They serve as prototypes for such behaviour in typical G2 solutions. In 
Chapter[3 the results of Chapters El and IHl are used to predict the past asymp- 
totic dynamics of G2 cosmologies. In particular, we discuss the existence of 
Mixmaster dynamics and the validity of the BKL conjecture II. 

Chapter |H1 is devoted to various aspects of numerical simulations of G2 
cosmologies: choosing the initial condition, devising the numerical solvers, 
and comparing the errors of the solvers. 

In Chapter |H1 we present the results from a variety of numerical simula- 



CHAPTER 1. INTRODUCTION 



11 



tions to illustrate various facets of G2 dynamics. In contrast to the previous 
numerical work on vacuum models, we investigate non-vacuum models, in 
which the source of the gravitational field is a radiation fluid and a cosmo- 
logical constant. The main focus is on the past asymptotic regime, with a 
secondary focus on the future asymptotic dynamics with a positive cosmo- 
logical constant and on the close-to-FL regime. 

In the concluding Chapter^] we summarize the new results in this thesis, 
and remark on future research. 

As regards background material, we assume that the reader is familiar 
with the orthonormal frame formalism, although a brief review with refer- 
ences is given at the beginning of Chapter |21 In addition we assume familiar- 
ity with some basic concepts in the theory of dynamical systems. We refer to 
WE (Chapter 4) for a summary and further references. For background ma- 
terial on the use of scale-invariant variables and dynamical systems methods 
in cosmology, we also refer to WE. 



Part I 

Dynamics of G2 and SH 
cosmologies 
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Chapter 2 



The orthonormal frame 
formalism 

In this chapter we introduce the orthonormal frame formahsm, which is based 
on a 1+3 decomposition of spacetime. Our goal is to apply the formalism to 
give a 1+1+2 decomposition of spacetime, which provides the basis for our 
analysis of G2 cosmologies. 

The orthonormal frame formalism was first introduced in cosmology by 
Ellis 1967 [201 and subsequently developed by Ellis & MacCallum 1969 [231 
and MacCallum 1973 [S21. More recently, it was extended by van Elst & 
Uggla 1997 [73J. Our standard references will be MacCallum 1973 ^2 and 
van Elst & Uggla 1997 

2.1 1+3 decomposition 

In the orthonormal frame approach one does not use the metric g directly 
(as done in the metric approach), but chooses at each point of the spacetime 
manifold (A^,g) a set of four linearly independent 1-forms {oj""} such that 
the line element can be locally expressed as ds^ = r]abU>"'UJ^, where rjab = 
diag(— 1, 1, 1, 1). The corresponding vector fields {e^} are then mutually 
orthogonal and of unit length - they form an orthonormal basis, with eg 
being timelike (and thus defining a timelike congruence). 

The gravitational field is described by the commutation functions 7^^^ of 
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the orthonormal frame, defined by 

[ea,efe]=7Vc- (2.1) 

The first step is to perform a 1+3 decomposition of the commutation func- 
tions as follows:^ 

[eo, e^,] = UaGo - [H5j + a^'^ - ej^iuj^ - ^]^)] e^g , (2.2) 
[e^, e^] = -2ec,/3^t^^eo + [e^fivn^'' + aJa^ - Q-a^d'] • (2.3) 

The variables in (|2.2p and (|2.3|) have physical or geometrical meanings, as 
follows. The variable H is the Hubble scalar, cTq,^ the rate of shear tensor, 
iioL the acceleration vector, and the rate of vorticity vector of the timelike 
congruence defined by eo, while ^q, is the angular velocity of the spatial 
frame {eo} with respect to a nonrotating frame (f^Q, = 0). The variables Hap 
and have no direct physical or geometrical meanings in general, but if the 
fundamental timelike vector field eo is hypersurface-orthogonal, then they 
directly determine the geometry of the spacelike 3-surfaces orthogonal to eo 
(see, for example WE, Section 1.6.3). We shall thus refer to n^^ and Oq, as 
the spatial curvature variables. Collectively, the variables above describe the 
gravitational field. We shall refer to them as the gravitational field variables, 
and denote them by the state vector 

-^grav (-^5 

The matter content of a cosmological model is described by the stress- 
energy tensor Tab, which is decomposed into irreducible parts with respect 
to eo in the following way (let eo = u below): 

Tab = fJ^UaUb + 2q(^aUb) + phab + T^ab, (2.5) 

where 

gaW" = 0, TTabU^ = 0, Ha" = 0, TTab = T^ba, 

and hab = Qab + UaUb is the projection tensor which locally projects into the 
3-space orthogonal to u. Since we are using an orthonormal frame, we have 
gab = Vab, = (1,0,0,0), and go = = vroa- The variables {^,p,qa,'Kaf3) 

^We follow the Via sign convention of van Elst & Uggla 1997 jZHl page 2676]. 



CHAPTER 2. THE ORTHONORMAL FRAME FORMALISM 



15 



have physical meanings: /i is the energy density, p is the (isotropic) pressure, 
Qa is the energy flux density and tTq/j is the anisotropic pressure (see, for 
example, van Elst & Uggla 1997 73j page 2677]). We shall refer to these 
variables as the matter variables, and denote them by the state vector 

Xmatter = {fJ', Qa, P, T^ap) ■ (2.6) 

The dynamics of the variables in ()2.4|) and ()2.6p is described by the EFEs, 
the Jacobi identities (using e^) and the contracted Bianchi identities: 

Gab + Agab = Tab , (2.7) 
eicYab] - Yelcl\b] = , (2.8) 

and VbTa'' = , (2.9) 

respectively. The evolution of p and TCap has to be specified by giving an 
equation of state for the matter content (e.g. perfect fluid). The variables iia 
and Qa correspond to the temporal and spatial gauge freedom respectively, 
and will be specified later. 

Throughout this thesis we shall assume that the fundamental congruence 
eo has zero vorticity {ua = 0), or equivalently that eo is orthogonal to a family 
of spacelike 3-surfaces. We now list the system of evolution and constraint 
equations (also see WE, page 33, but with the opposite sign convention for 
fla] see van Elst & Uggla 1997 [731 for the equations with vorticity). The 
angle bracket used in ()2.12j) and ()2.16|1 denotes traceless symmetrization - 
given Vai3, V^ap) is defined by 



(2.10) 
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The general orthonormal frame equations {ua = 0) " 

Einstein field equations: 

eoH = -H^ - laaf^a""^ + |(e„ + < - 2a«)«" - + 3p) + (2.11) 

+ {e^a + U{a + a{a)ui3) (2-12) 

= (Cg) = QH^ + ^R- a^f^a'^^ - 2^ - 2A (2.13) 

= (Cc)a = -2eaH + Bf^aJ - 3apaJ - ej^ripsa^' + , (2.14) 

where ^R and ^Sap are the isotropic and anisotropic spatial curvature, and 
are given by 

^R = 4ec.a" - 6a„a" - n^^n"'^ + | {n^''f (2.15) 

= e(c,a^) - (e^ - 2a^)ns^a^f3)^^ + 2n(«^n;3)^ - n^^ni^a/s) ■ (2.16) 

Jacobi identities: 

Goiriap) = -Hriap + 2a^^an/3)y - 2e'^^(«n^)^fi5 

- (e^ + u^){e^\aCTfs)5 - 5\a^p) + d^pW) (2.17) 
eo(aa) = {-HSj - a J - e„'^'^0^)a/3 

- i(e^ + U(,){2H5a^ - a J - e/^Q^) (2.18) 
= (CJ" = [£^^^^(6^ - a^) - n'^'^]u, (2.19) 
= (Cj)„ = ep{nj + e„^^a^) - 2a^n/ . (2.20) 

Contracted Bianchi identities: 

eo(/i) = -3//(/x + p) - ct^^tt"^ - (e« + 2< - 2a«)?'^ (2.21) 

eo(?a) = - a J - eJ'^^j)q/3 - {i^ + p)ua - ^aP 

- {ep + up- Zap)!: J + eJ-'upsT^^^ . (2.22) 



"The letter C denotes a constraint, and the subscript denotes which one: Gauss, Co- 
dacci, vorticity and Jacobi. 
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Perfect fluid 

This thesis focuses on cosmological models whose matter content includes 
a perfect fluid, which we now introduce. A perfect fluid has a 4-velocity u 
which is not necessarily aligned with the 4-velocity eg = u of our fundamental 
observers. The stress-energy tensor Tab is 

Tab = i^Uallb + P{gab + UaUb), UaU" = -1. (2.23) 

A perfect fluid with a linear barotropic equation of state specifies p: 

P = P{fi) = (7 - l)/i, (2.24) 

where 7 is a constant parameter. The range 

1 < 7 < 2 (2.25) 

is of particular physical interest, since it ensures that the perfect fluid satisfies 
the dominant and strong energy conditions and the causality requirement 
that the speed of sound should be less than that of light. The values 7 = 
1 and 7 = 1 correspond to pressure-free matter ("dust") and radiation, 
respectively. Since u and u are not necessarily aligned, we write 

ii!" = ^11" + (2.26) 

where 

u„ = 0, i;o = 0, r = (l-t;2)-i, v^ = v^v''. 
It follows from ()2.5|) . ()2.23|) - ()2.26p that {^,p,qa, 11^/3) are given by 

/i = -r^fi, (2.27) 

P = ^ qa = T^^Va, TTaf3 = 7^t^(a^^/3>, (2.28) 

where G+ = 1 + (7 — l)f ^. The basic variables that we use are /x and v". 
The vector f " is called the tilt of the fluid, and has three degrees of freedom. 
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Weyl tensor 

It is useful to keep track of the Weyl curvature tensor Cabcd, decomposed 
with respect to u = eo according to (see WE, page 19): 

and are given in terms of the gravitational and matter variables as follows 
(See WE, page 35): 

Ea/3 = HdafS - (o-a'^CTjlS - ^O'^SCr'^^ Sa/s) + ^Sa/S " ^7^0/3 (2.29) 

Haf3 = (e^ - a^)(rsia(^f3f^ - + n-^&a^^^ap + \n^'^aap . (2.30) 

Local coordinates 

For quantitative analyses, we need to express the differential operators 
eo and e^ as partial differential operators in terms of t and x*. ^ The 1+3 
decomposition for {eo,eQ,} is as follows: 

eo = N-^dt , e, = M^dt + ejdi . (2.31) 

The assumption that eo is hypersurface-orthogonal allows us to set 

M„ = 0. (2.32) 

It follows that the 3-surfaces^ to which eo is orthogonal are given by t = 
const.. See van Elst & Uggla 1997 73, Section 2.7] for details. The variables 

{N, ej) (2.33) 

are called the frame coefficients. When the commutators ()2.2|) and ()2.3|) act 
on t and x\ they provide the equations that relate the frame coefficients to 
the variables in (j2.4j) . 



^We reserve the indices i, j, fc, running from 1 to 3, for the local coordinate components 
with respect to x*. 

■^For brevity we shall refer to these 3-surfaces as "slices" . 
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Commutators acting on t and x*; 

eo(e,0 = i-H6j - - ej^^,)ep' (2.34) 
= iC^)a = N-'e^{N)-u^ (2.35) 

= (Ccom)'a/3 = 2(e[„ - a[„)e^f - ea(3sn"'^e^' . (2.36) 
The frame coefficients are directly related to the metric g through 

g00^_j^-2^ ^o.^Q^ ^^J- = 5-/5e„V. (2-37) 

This ties the variables in (|2.4p to the metric. 

All these equations give an autonomous system of first order PDEs involv- 
ing temporal and spatial derivatives, and a set of constraint equations which 
involve only spatial derivatives, for the variables ()2.4|) . ()2.6|) and ()2.33|) . The 
system is under-determined, however, since there are no evolution equations 
for N, iia, ^la, P and tTq^. A perfect fiuid matter source with linear barotropic 
equation of state gives ()2.28j) . while temporal and spatial gauges determine 
{N,Uo,) and fla respectively. 

The general orthonormal frame equations ()2.11|) - ()2.22|) describe the dy- 
namics of Go cosmologies. They are also ideal for describing the dynamics of 
SH cosmologies of Bianchi types VIII and IX, for which the t = const, slices 
are chosen to coincide with the group orbits (and the orthonormal frame 
is chosen to be group-invariant) , so that 

ea(Xgrav) = 0, ea(Xmatter) =0, ila = 0. (2.38) 

2.2 1 + 1+2 decomposition for G2 cosmologies 

By a G2 cosmology, we mean a cosmological model which admits an Abelian 
group G2 of isometrics whose orbits are spacelike 2-surfaces. We derive the 
evolution and constraint equations for G2 cosmologies by specializing the 
general orthonormal frame equations in the previous section. We shall pre- 
serve the vorticity-free eg while taking full advantage of the G2 structure to 
simplify the equations, namely to make the variables depend on t and x only. 

Let ^1, ^2 be two independent Killing vector fields (KVFs) of the Abelian 
group of a G2 cosmology. First, we choose a group-invariant orthonormal 
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frame, i.e. 

[e„,ei] = 0= [e„,e2]. (2.39) 

Secondly, we also adapt the orthonormal frame {6^} so that 62 and are 
tangent to the G2 orbits. The existence of such a frame is established in the 
proof of Theorem 3.1 in Wainwright 1979 j73] (see Hewitt 1989 [34, page 120] 
for a more general theorem). For a group-invariant orbit-aligned frame it 
follows (Wainwright 1979 (TS]) that 

eA{l%c) = 0, (2.40) 
7%A = 7V = 0, 7V = 7^1^ = 0, 7023 = 7123 = 0. (2.41) 

Recall that Latin indices run from to 3; Greek indices from 1 to 3; capital 
Latin indices run from 2 to 3. Writing (|2.41|) in terms of the irreducible 
variables in ()2.4p gives 

iiA = 0, Qa = -^AB<y^^, = 0, aA = 0, riia = 0. (2.42) 

Note that Qa is now expressed in terms of ctia, and eAs is the anti-symmetric 
permutation tensor (e23 = —632 = 1). 

For convenience we introduce the trace-free quantity 

nAB = nAB - \nc^5AB, (2.43) 

and its dual 

*nAB = eAcnB^, (2.44) 

and similarly for cab, ^Sab and tt ab- The resulting system involves the 
following variables: 

Xgrav = {H, (Til, CTlA, ^AB, , flAB, Ml, ^^l), (2-45) 

Xmattcr = (yU, 71^11, T^IA, Ab) ■ (2.46) 

By performing a 1-1-2 spatial decomposition of the general orthonormal 
frame equations in the previous section, we obtain the following equations 
for G2 cosmologies.^ These equations arise in a fairly obvious way from the 
general orthonormal frame equations ()2.1H) - ()2.22|) . with VLa replaced using 



^ These equations generalize the system of equations in WE, pages 48-49, which are 
vahd for OT G2 cosmologies. 
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The 1+1+2 equations for G2 cosmologies 

Einstein field equations: 

eoH = -H^ - \o^pG'^^ + |(ei + ii^ - 2a^)u^ + 3p) + |A (2.47) 

eo(f^ii) = -3i/aii - Ig^a^^^ + |(ei + mi + ax)ux - '^Syv + tth (2.48) 
eo(f^iA) = -3(i? - \axx)(yxA - {^ab - ^leAs^^ + ttia (2.49) 

eo(5-AB) = -SHgab + 2criA(TiB - GicCT^'^Sab - 2fli*aAB - U\flAB 

-'^AB+^AB (2.50) 

= (Cg) = QH"^ + ^R- a^pa''^ - 2/x - 2A (2.51) 

= (Cc)i = -2ei{H - ian) - SaKm - *aABn^'' + Qi (2.52) 

= {Cc)a = (ei - 3ai)(7iA - {^uab + |nc^eAB)(7^^ + qa , (2.53) 

where 

fTa/3(7"'^ = i((Jll)' + 2(71^^1^ + aAB^^^ (2.54) 

^i? = 4eiai - 6a? - UABn^^ (2.55) 
%i = feiai - fnAsn^^, %a = (2.56) 
= (ei - 2ai)*fiAB + nc^UAB ■ (2.57) 

Jacobi identities: 

eoiriA^) = (-/^ - (Jii)nA^ + 2a^^nAB - 2(ei + -^1)^1 (2.58) 

eoifiAB) = (--H" - crii)nAB + uc^uab - 2Qi*fiAB + (ei + iixfoAB (2.59) 

eo(ai) = {-H - (7ii)ai - (ei + u^) {H - ifjn) . (2.60) 

Contracted Bianchi identities: 

eoin) = -3H{fjL+p) - fcriiTrii - 2aiAn^^ - ^abtt^^ 

- (ei + 2ui - 2ai)qi (2.61) 
eo(?i) = (-4i? - (7ii)gi - 2(7iAg^ - eip - {ii + p)ui 

- (ei + -fii - 3ai)7rii - *nABT^^^ (2.62) 
eo(?A) = (-4i7 + |(7ii)gA - (5-ab + f^ie^s)?^ 

- (ei + Til - 3ai)7riA + Cuab + hnc^ eAB)n^^ ■ (2.63) 
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Perfect fluid 

The 1+2 spatial decomposition of ()2.28|) gives 

(7-l)(l-i;2) + |7t;2 
P = ^ Qi = -p^n, Qa = -t^va, 

(jr_|_ ^+ ^ + 

(_T-[- ^ + 

ttab = yrk^AVB - Ivcv'^6ab) , (2-64) 



where 



= 1 + (7 - = v^v» = yj + v^v^, (2.65) 



VKey/ tensor 

The 1+2 spatial decomposition of ()2.29|) - ()2.30|) gives 

En = i/an - l^n' - + la^B^'^'' + '^n - (2.66) 

EiA = HaiA - {\(JiiOiA + o^^gab) - \t^ia (2.67) 

Eab = Hoab - (cTiyiais - \aicCi^^5AB - dnaAB) + oab - \t^ab (2.68) 

i/ii = UAB^'^'' (2.69) 

iJiA = -\{d^ - ai)aii?eA'' - l^'^'^AB - i^c^c^iA (2.70) 

^AB = -(^1 - ai) *S^AB - nc'^CTAB + fo^ii^-AB • (2.71) 

Local coordinates 

We choose the x-coordinate to be constant on the G2 orbits. The 1+1+2 
coordinate decomposition of the frame coefficients (|2.33|) is as follows: 

eo = N-'^dt, ei = ci^d^ + = ca^O^ + ca^O^i } (2.72) 

That X is constant on the G2 orbits means ga^x) = 0, which implies 

ca' = 0, (2.73) 



^We reserve the indices /, J, K, running from 2 to 3, for the local coordinate components 
with respect to y and z. 
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It follows from (jOOl), ff^^ and (jT7H|l that 

9,.(Xgrav) = 0, 9,.(X^atter) = 0, (2.74) 

which implies that in the evolution and constraint equations the operator ei 
is effectively 

ei = ei^^, . (2.75) 

In conjunction with ()2.42|) . the commutators ()2.2|) and ()2.3p acting on t, x 
and give the following equations for the remaining frame coefficients 

(iV,ei\e/,e/). (2.76) 

Commutators acting on t, x and xh 

= (a)i = iV-iei(iV)-iii (2.77) 

eo(eii) = (-iJ-(Tn)eii (2.78) 

eo(e/) = {-H - (Tii)e/ - 2a^^eA^ (2.79) 

eo{e/) = -He/ - aA^'cB' - nieA^'es' (2.80) 

= (Ccom)'iA = (ei - a,)e/ - (*n^^ + inc^e^^)eB'. (2.81) 

Note that ei^ and ca^ decouple from the system ()2.47|) - ()2.63|) . 

Equations dTTTD - dZTK^ . (ETTfl and (IT7HD give a system of first order 
PDEs in t and x for the variables in (j2.4j) - (j2.6|) and the pair (A^, ei^). The 
system is under-determined, however, since there are no evolution equations 
for N, ill, Qi, p and Hap- A perfect fluid matter source with linear barotropic 
equation of state gives ()2.64|) . while temporal and spatial gauges determine 
{N,ui) and fli respectively. These equations are used to analyze the dynam- 
ics of G2 cosmologies and of SH cosmologies whose G^, group has an Abelian 
G2 subgroup. 



2.3 Scale- invariant variables 

The variables in (|2.4j) . (|2.6|) and (|2.33|) are scale-dependent and dimensional, 
and are unsuitable for describing the asymptotic behaviour of cosmologi- 
cal models near the initial singularity, since they typically diverge. It is 
thus essential to introduce scale-invariant (dimensionless) variables, which 
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one hopes will be bounded as the initial singularity is approached. As a 
motivating example, consider an FL model, in which the matter density 
/i — >• cx) at the initial singularity. One wants to introduce a correspond- 
ing scale- invariant, dimensionless density parameter Q = /i/(3if^) such that 
Q ^ 1.^ Another example is the Kasner models, in which the shear cTq,^ sat- 
isfies cr^ = ^<7ai3<7°'^ ^ oo at the initial singularity. One wants to introduce a 
corresponding shear parameter T,a/3 = o-a/s/H such that = o"^/(3iJ^) — > 1. 
These examples motivate the use of the Hubble scalar H as the normalizing 
factor. 



1+3 Hubble-normalized variables 

We define the Hubble-normalized gravitational and matter variables respec- 
tively as follows: 

{T,af3,Ua, Ra, Na/3, Aa) = {(Ja/3 , Ua, ^a, ^a/S , tta) / H (2.82) 

(fi, Q„, P, H,^, ^a) = (/i, g„, p, TT^fs, A)/i3H^) . (2.83) 

The differential operators ()2.31|) and the frame coefficients ()2.33|) are also 
redefined to be scale-invariant: 

^0 = -^eo, da, = -^e„, (2.84) 

X = NH, E^' = ^, (2.85) 

since the frame variables have the same dimension as H. In order to trans- 
form the evolution and constraint equations, it is necessary to introduce the 
deceleration parameter q and the Hubble gradient according to 

doH 

q + l = -^ (2.86) 

(2.87) 



H 

In order to transform the equations, it is convenient to use the following 
identities. For the gravitational field variables Xgrav in ()2.4|) . which have the 



^Although it is sufRcient to require that fl tends to a non-zero constant, it is customary 
to normalize /i such that 17 — s- 1. 
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same dimension as H, we use the identities 



^,eoX,,,. = [^0 - (g + 1)] [^j (2.88) 
-^e.Xg,,. = [d^ - (^^^ , (2.89) 

which follow from ()2.82|) . ()2.84p - ()2.87p . Similarly, for the matter variables 
Xmatter in (EISl), which have the dimension [H"^], we use the identities 

^eoX^a^er = [do - 2(g + 1)] (^^) (2-90) 
^ re„X^,tter = [d. - 2r J ( • (2.91) 



f^^;^3 « matter a ctj i ^^^^ 

Then equation ()2.11|) for d^H decouples, and through ()2.86|) gives the fol- 
lowing expression for q: 

q = + \{n + 3P) - - \{d^ -r^ + Ua~ 2A^)U", (2.92) 

i.e. q is expressed in terms of the Hubble-normalized variables ()2.82|) and 
Equations (j2.28p in Hubble-normalized form are 

(7-l)(l-^^) + >^ ^ n n ^9 0-^^ 

" = 7^ ^i, Qa = -^Va, ilaf3 = -^^{^Vp) ■ (2.93) 

The spatial curvature variables % and '^8^13 are normalized according to 

n,.-^, 5„,^||. (2.94) 

The Weyl curvature variables Eap and iJo/j are normalized according to 

{£ap, Uap) = {Eap, H^p)/{'iH^). (2.95) 
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The commutators ()2.2j) - ()2.H|) in Hubble-normalized form are 

[^0, d^] = {U^ - r^)d^ - {-qda" + + eJ^R^) dp , (2.96) 
[5a, dp] = [e„^.iV^^ + (A, + r,)5/ - {Ap + r^)5„^] 8^ . (2.97) 

The 1+3 Hubble- normahzed variables and equations are well-suited for 
analyzing Gq cosmologies. They are also well-suited for SH cosmologies of 
types VIII and IX, in which (|2.38|) implies 

a,(Xg,av/i^) = 0, a„(X^atter/(3i/')) = 0, = 0, r„ = 0. (2.98) 

It then follows from and that we can set 

A/" = 1, so that do = dt . (2.99) 

1+1+2 Hubble-normalized variables for G2- compatible 
SH cosmologies 

For G'2-compatible SH cosmologies, we define the Hubble-normalized 1+1+2 
decomposed gravitational and matter variables as follows: 

{T,ii, T,iA,t,AB, R, Nc^ , Nab, A) = {<yii,aiA,aABMi,nc^ ,nAB,ai)/H 

(2.100) 

{Q,Qi,QA,P,'n.ii,'n.iA,'n.AB,^A) = (//, gi,gA,P, 7rii,7riA,vrAB,A)/(3iI^). 

(2.101) 

Equivalently, one can view this as the 1+1+2 decomposed version of the 
1+3 Hubble-normalized variables in (EH -(EH- Equations 
are similarly decomposed. 

1+1+2 /3-normalized variables for G2 cosmologies 

For G2 cosmologies, the area expansion of the G2 orbits is described by 

eAB = H5AB + (TAB. (2.102) 

The average expansion rate of the G2 orbits is described by 

/3 = iec^ = if-ian. (2.103) 
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It turns out that there are some advantages in using (3 instead of H as the 
normahzing factor (see van Elst et al. 2002 f73]). We define the /3-normahzed 
variables as follows: 

{T.n,T.iA,tAB,U , R, Nc^ , Nab, A) = {(Tu, criA,^AB, ui, ^i, nc'^ ,nAB,ai)/ P 

(2.104) 

(n,Qi,QA,^, nii,niA,nAB,^A) = (/i,gi,gA,P,7rii,7ri^,7rAB,A)/(3/5^). 

(2.105) 

The differential operators are likewise normalized as follows: 

9o = ^eo, di = ^eu (2.106) 

Af = N/3, ^i' = y, (2.107) 

We define the analogous deceleration parameter q, and the /5 gradient r 
according to 

«+l = -f. r = -f. (2.108) 

The equations are similarly transformed using the following identities: 

i^eoX,,,. = [^0 - (g + 1)] (2. 109) 

le,X,,.. = [5i - r] (2.110) 
^eoX^,,,er = [^0 - 2(g + 1)] (^^) (2.111) 

1 /V 

„ 1 / -^matter 



3^3eiX^.,er = [d, - 2r] ^^^J . (2.112) 

Consider the evolution equation for P, in which eiOi is replaced using (Cq): 

eoP = + \aiAcy^^ - \{&ab^^^ + uabu^^) + \al - aiiii 

-\{p + 'nii) + \A. (2.113) 
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This defines the corresponding area deceleration parameter q: 

g = i - iSi^S^^ + + NabN^"") - \A^ + AU 

+ |(P + nn)-|^^A. (2.114) 

The evolution equation for 13 decouples from the main system. The (Cc)i 
constraint ()2.52p can be solved for r: 

r = |(ASn + *tABN^'' - Qi). (2.115) 
The spatial curvature variables ^Su and ^Sab are normalized according 

to 

'^R "^Sii ^Sab 

The Weyl curvature variables are normalized according to 

{£ii,£iAi£abi'Hii,1-Lia^'Hab) = {En, Eia, Eab^ Hu, Hia, Hab) / {'iP'^)- 

(2.117) 

The commutator [^oj^i] is given by 

[^c^i] = {U- r)do + (g - |Sn)5i - 2S/5^ . (2.118) 



Chapter 3 

The Go and SH hierarchies 



In this chapter we give a classification of G2 cosmologies, and of the SH 
cosmologies that arise as limiting cases, using the action of the G2 group of 
isometrics, and the tilt degrees of freedom of the perfect fluid. Our goal is to 
provide a framework for describing the remarkably diverse range of dynamical 
behaviour that is possible within the class of G2 cosmologies, only a small 
part of which has appeared in the literature to date. 

3.1 The G2 hierarchy 

It is useful to classify G2 cosmologies according to the action of the G2 
group, as determined by the properties of the KVFs associated with the G2 
group. There are four restrictions relating to the KVFs, which determine four 
special classes of G2 cosmologies. If none of these restrictions is satisfied we 
shall refer to the G2 cosmology as generic. We now give the above-mentioned 
restrictions, and then characterize the subclasses using the gravitational field 
variables 

{l)Hypersurf ace- orthogonal 

This case is specified by requiring that the G2 admits a hypersurface- 
orthogonal KVF ^, 

i[a±ic] = 0. (3.1) 

The following proposition provides a convenient way to check for the existence 
of a hypersurf ace-orthogonal KVF. 
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Proposition 3.1. Consider a G2 cosmology presented in a group-invariant 
orbit-aligned frame with 62 and tangent to the G2 orbits. Then the G2 
admits a hypersurf ace- orthogonal KVF if and only if the gravitational field 
variables \2.4^ satisfy 

<^AB^ = U, ^Aficr (J = U, crAB<^ cr = (J. ['J-^j 
In addition, the matter variables ^^^J^ satisfy 

e^'^QAcyiB = 0, e^'^'K.AcriB = 0, *^AB(y^^(y''' = 0. ' (3.3) 
Proof. See Appendix |X| □ 



(2) Orthogonally transitive 

This case is specified by requiring that the G2 acts orthogonally transi- 
tively, 

i[a;bicrid] = 0, ri^a^hricid] = 0, (3.4) 

where ^ and rj are two linearly independent KVFs (Carter 1973 fSl page 
160]). In other words, the 2-spaces orthogonal to the orbits of the G2 are 
surface-forming. This is equivalent to (see WE, page 46) 

a,A = 0. (3.5) 

In this class, equations ()2.49|) and (j2.53|) imply that the matter variables 
()2.6j) are restricted as follows: 

qa = 0, TiiA = 0. (3.6) 

(3) Diagonal 

This case is specified by requiring that the G2 admits two hypersurface- 
orthogonal KVFs. This class is the intersection of the first two classes, and 



^If fiAB and a^^ are non-zero, then any two of the conditions in (|8.2(l imply the third. 
Furthermore, for SH cosmologies, we also have nc'~^ — 0. 
^It also follows from H3.2|l and that 



^Asa^"" = 0, T^Asn^^ = 0. 
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is thus characterized by 

(TiA = 0, aABn^"" = 0, (3.7) 

as follows from ()3.2j) and ()3.5|1 . In this class the matter variables ()2.fij) are 
restricted as follows: 

qA = 0, 7r,A = 0, *nAB^^^ = 0. (3.8) 

The name "diagonal" reflects the fact that the line element can be written 
in diagonal form in this class. 

(4) Plane symmetric or locally rotationally symmetric (LRS) 

In the diagonal G2 class, if there exists in addition a one-parameter 
isotropy group, then the model is locally rotationally symmetric (WE, page 
43). This class is characterized by 

o-iA = = aAB = riAB , (3.9) 

since otherwise there would exist a preferred direction in the 23-space. In 
this class the matter variables ()2.6|) are restricted as follows: 

Qa = 0, TTiA = 0, ttab = 0. (3.10) 

Figure summarizes the classification of G2 cosmologies using the action 
of the G2 group and the existence of an LRS group (see WE, page 44). On 
account of ()3.2p and ()3.5p we can say that a G2 cosmology is generic if and 
only if 

a^Aa'^ + 0, [dABU^^'f + {uABO^^o^^'f + {:aABO^^a'''f ^ 0. (3.11) 




Figure 3.1: Classification of G2 cosmologies. 



Assuming a perfect fiuid matter content, we now show that G2 cosmolo- 
gies can be further classified using the tilt degrees of freedom of the per- 
fect fiuid, as described by the tilt vector Va (see Section EH}. From ()2.49|) . 
fl2.63|) . and using ()2.53|) and ()2.64|) . one finds that the evolution equation for 
e^^g^CTiB is homogeneous: 

eo(e^^9AO"iB) = {-IH + 2aii-(iiVi-ViUi + Qviai-Viei){e'^^ qAOiB)- (3.12) 

Recalling from (j2.64|) that = iG^^^va, equation (j3.12|) implies that 
e^^^VA^yiB = defines an invariant set. Likewise, the evolution equation 
for qA is homogeneous, on account of ()2.64|) . implying that f a = forms 
an invariant set. These considerations lead to the hierarchy presented in 
Table 



■^By inspection of equations (|2.62() - (|2.63() and using H2.64|l . one finds that cases where 
fi = 0, f A 7^ are not invariant. 
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Table 3.1: Perfect fluid hierarchy in the presence of a G2 group. 



Restriction 


Description 


tilt degrees of freedom 


e^^^A^TiB ^ 


va not parallel to ctia 


3 


e^'^VACTiB = 0, va^O 


va parallel to aiA 


2 


va = 0, vi^O 


Va orthogonal to G2 orbits 


1 


Va = 


non-tilted fluid 






Figure|ni21summarizes the classification of G2 cosmologies using the action 
of the G2 group and the tilt degrees of freedom. Note that the matter restric- 
tion ()3.3j) of the class of hypersurface-orthogonal G2 implies e^^VAcriB = 0, 
and restriction ()3.(j|l of the class of OT G2 implies va = 0. The class of 
hypersurface-orthogonal G2 also contains va = as an invariant subset. 

3.2 G2-compatible SH cosmologies 

Perfect fluid SH cosmologies of Bianchi type I-VII admit an Abelian G2 sub- 
group (this can be inferred from Table 11.2 in MacCallum 1979 The 
properties of det uab and ai for G'2-compatible SH cosmologies are summa- 
rized in Table 13.21 These cosmologies can be further classified according to 
the G2 action (see Figure IXT|) and the tilt degrees of freedom, as in Table ITTl ^ 



Table 3.2: G'2^compatible SH cosmologies. 





ai = 


ai ^ 


det Uab > 


VIIo 


VII, 


det Uab < 


VIo 


VI, 


det Uab = 


II 


IV 


Uab = 


I 


V 



^The parameter h in Table IXTl will be introduced in equation (|3.19|l . 
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generic G2 action 



with 








with 









with va — ^ 



one HO KVF 



with e vaOxb — 



with va — ^ 



OT G2 {VA = 0) 




diagonal G2 {va — 0) 



Figure 3.2: Action of the G^ group and the tilt degrees of freedom. 
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SH cosmologies with generic G2 



Consider the SH cosmologies with a generic G2 subgroup. We now show 
that the number of tilt degrees of freedom imposes strong restrictions on the 
possible Bianchi type, and that conversely, for a given Bianchi type, the tilt 
degrees of freedom are restricted. The results are summarized in Figure 13.31 
The number on the right of each box in Figure 13.31 indicates the tilt degrees 
of freedom, as defined in Table 13.11 



det uab 0, af = h det uab det uab = 



VI, VII 



VI 



VI* 1 



VP 1 



IV 



Uab = 



IV 



V 



II 



Figure 3.3: SH cosmologies with generic G2 action. The numbers 0,1,2,3 
indicate the tilt degrees of freedom. 



We now justify the results in Figure 13.31 Bianchi V cosmologies with 
generic G2 action necessarily have aiA parallel to qa, as follows from (|2.53p . 
Note that Bianchi VII cosmologies with generic G2 action necessarily have 
the most general tilt, as shown in the following proposition. 



Proposition 3.2. Consider a G 2- compatible SH cosmology presented in a 
group-invariant orbit-aligned frame with 62 and tangent to the G2 orbits, 
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and ei tangent to the G3 orbits. If 

det uab > and gxao^^ 7^ 0, 

then 

e^^'vAcriB ^ 0. 

Proof. Suppose e^^VAC^iB = 0. Contract (Cc)a with e^^^aiB to obtain 

nc%Acr^^ = 2hABCT^^cr^^- (3.13) 
Assuming aiAO'^^ 7^ 0, represent aiA using polar coordinates as follows: 

(JlA = (7±KiA, 

where K12 = cos 0, K13 = sin 0. Then it follows from ()3.13|) that 

nc^ = 2fiABK^^K^'^ ■ (3.14) 

The determinant of hab-, given by 

det(n^B) = {{nc'^f - InABn^"" , (3.15) 

is then expanded and simplified using ()3.14|1 and trigonometric identities to 
yield 

det(n^B) = -iCuABK^^K^^'f < 0. (3.16) 
where *hAB is defined by ()2.44|) . □ 

The h parameter for Bianchi Ylh and VII/^ cosmologies 

We now define the h parameter for Bianchi Vlh and Vllh cosmologies. Note 
that the evolution equation for det(nAB) is given by 

eo(detnAB) = 2(-/7 - crn) det n^B, (3.17) 

as follows from ()3.15|1 and ()2.58j) - ()2.59j) . Comparing with the evolution equa- 
tion for ai: 

eo{ai) = {-H - au)ai, (3.18) 
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it follows that if det(n^B) 7^ 0, there exists a constant h such that 

al = hdet{nAB)- (3.19) 

Bianchi VII^ cosmologies are characterized by 

det(nA_B) > 0, cti 7^ 0, 

requiring that h satisfies h > 0, while Bianchi Vlh cosmologies are charac- 
terized by 

det(?2A_B) < 0, ai 7^ 0, 
requiring that h satisfies h < 0. 

The exceptional Bianchi VI% cosmologies 

9 

We now explain the origin of the exceptional Bianchi Vl% cosmologies. For 

9 _____ 

Bianchi Vlh models that also satisfy e^^VACis = 0, we can use ()3.16p to 
write 

ai = ±v^ *hABK^'^K^^- (3.20) 
Notice that from {Cc)a(t^^ yields 

= -3[ai + \*nABK'''K'^] + (3.21) 

for cr_L 7^ 0. Thus, choosing ai + ^*hABK^^K^^ = 0, which results in h = — |, 
also forces QaK^^ = 0. In conjunction with e^^VAO'iB = and cr_L 7^ 0, this 
yields 

va = 0. (3.22) 

This class of models is referred to as the class of exceptional Bianchi Vl% cos- 

9 

mologies. In other words, the exceptional Bianchi VI% cosmologies are char- 

9 

acterized by a non-OT G2 action {a± 7^ 0) and a fluid congruence that is or- 
thogonal to the G2 orbits {va = 0). Note that choosing ai — Y'it-abK^^K^^ = 
also results in h = — |, but does not permit qaK^^ = 0. This class consists 
of the ordinary Bianchi Vl_i cosmologies. 
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Bianchi IV and II cosmologies 

Note that it follows from and (IT^ that 

det uab = 0, va = 

imphes 

Oi = 0, 



(3.23) 



(3.24) 



i.e. there is no Bianchi IV with generic G2 and one tilt degree of freedom. 
Lastly, we show that 



det riAB = 0, ai = 



Va = 0, 



(3.25) 



i.e. Bianchi II cosmologies with generic G2 have one tilt degree of freedom 
only. To see this, choose the spatial gauge ()4.65|) . Then det hab = and 
fl2.53|) imply that A^x = and va = 0. But f a = is rotation-independent, 
so it is not an artifact of the spatial gauge choice. 

Bianchi cosmologies with a generic G2 subgroup are interesting because 
all of them except Bianchi V exhibit the Mixmaster dynamics, as we shall 
show in Section 14.11 

Bianchi cosmologies with non- generic G2 

For completeness, Bianchi cosmologies with special G2 action are classified 
in Figures I3.4H3.51 



IV, VI, VII 



II 



v y_ 



IV, VI, VII 




n 


- 



Figure 3.4: SH cosmologies with OT G2 action. The numbers 0,1 indicate 
the tilt degrees of freedom. 
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-) I I I b) 

^ v,vi , : ^ 



V 



V, VI 




1 


• 


V (LRS) 


1 














I, V, VI 





;:::> 


I (LRS) 






Figure 3.5: a) SH cosmologies whose G2 admits one HO KVF. b) SH cos- 
mologies with diagonal G2 action. The numbers 0,1,2 indicate the tilt degrees 
of freedom. 



Chapter 4 
SH dynamics 



Spatially homogeneous cosmologies are relevant in the study of G2 cosmolo- 
gies because their dynamics serves as the background dynamics of the inho- 
mogeneous cosmologies in regimes in which the spatial derivatives are small 
compared with the time derivatives. In this chapter, we discuss certain as- 
pects of the asymptotic dynamics of SH cosmologies, in particular the occur- 
rence of Mixmaster dynamics. We introduce a unified system of evolution 
equations for the G2 compatible SH cosmologies (models of Bianchi types 
I- VII), which enables us to predict in which classes of models Mixmaster 
dynamics will occur. We also mention some of the unsolved problems con- 
cerning the dynamics of SH cosmologies, which can be addressed using these 
evolution equations. 



4.1 The standard Mixmaster dynamics 

The standard Mixmaster dynamics arises in vacuum SH cosmologies of Bianchi 
types VIII and IX. The justification for considering vacuum cosmologies is 
that for typical SH cosmologies with a non-tilted perfect fluid, the matter 
density is dynamically insignificant in the sense that 

lim ^] = 0, (4.1) 

where t is an appropriately defined time. (See Ringstrom 2001 [H^ for a 
proof of this result for the class of non-tilted Bianchi type IX cosmologies.) 
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The evolution equations 

The evolution equations for vacuum Bianchi VIII and IX cosmologies can be 
derived from the system ()2.11|) - ()2.22|) by choosing a Gs-invariant orthonor- 
mal frame and introducing Hubble-normalized variables as in Section 12. HI 

Recall from Table 11.11 that Aa = for models of Bianchi types VIII and 
IX. We use a time-dependent rotation of the spatial frame vectors {ea} to 
diagonalize Nap: 

Na/3 = diag{Nn,N22,N33). (4.2) 

The (Cc)i constraint ()2.14j) implies that the shear is also diagonal. It then 
follows from (j2.17|) that the spatial frame vectors are Fermi-propagated, i.e. 
Ra = 0. As in WE, Section 6.1.1, we introduce the notation 

S+ = |(S22 + S33), S_ = ^(S22 - S33). (4.3) 

The state vector ()2.82|) reduces to 

X=(S+,S_,iVn,iV22,iV33). (4.4) 

The evolution equations for X are obtained by taking the appropriate 
linear combinations of (j2.12|) and by specializing (j2.17|) : 

= {q- 2)S+ - I [(iV22 - N,,f - Nn{2Nn - N22 - N,^)] (4.5) 
= {q- 2)S_ - ^(iV33 - iV22)(iVn - iV22 - iVgs) (4.6) 

iV{, = (g-4S+)iVn (4.7) 
N^2 = (9 + 2S+ + 2v^S_)iV22 (4.8) 
iV^3 = {q + 2S+ - 2v^S_)iV33 , (4.9) 

where ' denotes the time derivative, and 

g = 2S^ S2 = S^ + S?. . (4.10) 

The only remaining constraint is the Gauss constraint ()2.13|1 . which reads 

= (CG) = l-fi;^-S^ (4.11) 
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where 

Qk = ^ [N'u + Ni, + iV^g - 2(iVniV22 + + N,,N^,)] . (4.12) 

The Kasner circle 

The Mixmaster dynamics is generated by the stabihty properties of the Kas- 
ner solutions. These vacuum solutions appear as a circle of equilibrium points 
of the evolution equations ()4.5|1 - ()4.9|) . given by 

A^„„ = 0, S^ + S'=l. (4.13) 

It is customary to refer to this circle as the Kasner circle, denoted /C. It 
follows from ()4.1()j) and ()4.13|1 that the Kasner solutions satisfy 

q = 2. (4.14) 

We note that the metric for the Kasner solutions has the simple form: 

e-^ds^ = -dT^ + T^Pirfx^ + T^'P^dy^ + T^P^dz^, (4.15) 

where the constants pa satisfy 

Pi +P2 +P3 = 1 = +P2 +P3- (4-16) 

The link between the pa and the Hubble-normalized shear is (see WE, equa- 
tion (6.17)) 

Pi = |(l-2S+), p2 = |(l + S+ + v^S_), p3 = |(l + S+-y3S^). (4.17) 

A typical Kasner solution is in fact represented six times on the Kasner 
circle. Figure l^?Tk illustrates this symmetry of the Kasner circle. The grey 
dots connected by dashed lines are six different representations of a typical 
Kasner solution. The exceptions are the so-called Taub Kasner solutions 
{pa = (1,0,0) and cycle), labelled on the Kasner circle in Figure HH^, 
and the LRS Kasner solutions {pa = (— |, |, |) and cycle), labelled Qa- See 
WE, Sections 6.2.2 and 9.1.1 for details. 
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Figure 4.1: a) The Kasner circle and its six-fold symmetry, b) The triggers 
^11, N22 and and their unstable Kasner arcs. 

Kasner instability and Mixmaster dynamics 

In order to proceed with the analysis we restrict our consideration to so- 
lutions that are close to a Kasner solution at some time, i.e. to solutions 
whose orbit in the Hubble-normalized state space enters an arbitrarily small 
neighbourhood of the Kasner circle /C. Ringstrom 2001 [HII proved that all 
Bianchi IX solutions except a set of measure zero have this property. 

Consider linearizing (j4.7j) - (j4.9|) on the Kasner circle (j4.13|) . Using (j4.14p . 
we obtain 

iV(, = 2(1 - 2S+)iVn (4.18) 
Ar^2 = 2(1 + S+ + y3S_)Ar22 (4.19) 
iV^3 = 2(l + S+-v^S_)iV33 , (4.20) 

which imply that the N^a are unstable into the past on the following arcs: 

Nil on arc(T2T3), N22 on arc(T3Ti), A^33 on arc(TiT2). 

These arcs are non-overlapping and cover all the Kasner equilibrium points 
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except the Taub Kasner points (see Figure HHJo). As a result, each of these 
Kasner equihbrium points is unstable into the past. In the SH state space, 
the unstable manifold of a Kasner equilibrium point is an orbit lying on one 
of the ellipsoids 

1 = n^L + + , « = 1, 2 or 3, (4.21) 

as follows from the Gauss constraint ()4.11|) . In fact, the orbits represent 
the Taub vacuum Bianchi II solutions, and appear as straight lines when 
projected onto the (S_|_, S_) plane, as illustrated in Figure lOl We shall refer 
to these orbits as the curvature transition sets, denoted T/v- Each transition 
is a result of the growth (into the past) of one of Naa- Therefore we call Naa 
the "triggers". Because Naa represent the spatial curvature, we call them 
the "curvature triggers" . 



^3 




B2 

Figure 4.2: The curvature transition set and the arc of instability (into the 
past) of the trigger A^n. The transition sets for and N33 are obtained by 
symmetry, using the points B2 and B3 respectively. 
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The evolution of a Bianchi VIII or IX cosmology is represented by an orbit 
in the Hubble-normalized state space. As t ^ — oo, a typical orbit enters a 
neighbourhood of the Kasner circle, but cannot be asymptotic to a Kasner 
equilibrium point in general. Instead it shadows the unstable manifolds of 
the Kasner equilibrium points, going from one Kasner equilibrium point to 
another. In this way, the three families of curvature transition sets sustain 
an endless sequence of transitions, called Mixmaster dynamics, in recognition 
of the work of Misner 1969 [Sl].^ For a more detailed discussion, see Ma & 
Wainwright [Sl^ and Section 6.4 of WE. The essential point is that Mixmaster 
dynamics will occur if for each typical Kasner equilibrium point P there is a 
trigger variable that renders it unstable into the past, and with the resulting 
unstable manifold joining P to another Kasner equilibrium point. 

To summarize, the occurrence of Mixmaster dynamics and the identi- 
fication of the past attractor is established by proving three results (see 
Ringstrom 2001 [Ulj): 

i) lim = 0, lim VL^ = 0. (4.22) 

ii) A typical orbit enters an arbitrarily small neighbourhood 
of the Kasner circle, and 

iii) lim {NuN22,N22N33,N3^Nn) =0. (4.23) 

f— >— oo 

Mixmaster dynamics is initiated by Property ii), which we shall refer to as 
the Kasner Attractivity Conjecture. ^ It appears to play a fundamental role 
in the dynamics of SH cosmologies and more generally, in inhomogeneous 
cosmologies, although it has not been proved in more general situations. 
Property iii) states that the curvature transitions cannot occur simultane- 
ously. 

For Bianchi VIII and IX cosmologies, the triggers are the eigenvalues of 
Na/3. If one of the triggers is zero, e.g. A''ii = 0, then there is a stable arc on 
K, which is the past attractor. For G'2-compatible SH cosmologies (Bianchi 
I- VII), at least one of the eigenvalues of Na/s is zero (see Table HHj). Can 
Mixmaster dynamics occur in these cosmologies? At first sight, it does not 



^Misner discovered this behaviour by using the Hamiltonian formulation of the EFEs. 
^This is stated formally later as Conjecture 14. II 
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seem so. But it has been shown that in fact it can occur, although in general 
it requires the presence of a perfect fluid with non-zero tilt, as we shall show 
in Section 031 But first we need to derive the Hubble-normalized evolution 
equations for G'2-compatible SH cosmologies, which we do in Section 14.21 

4.2 The G'2-compatible SH evolution equations 

The 1+1+2 Hubble-normalized variables ()2.100|) - ()2.10H) were introduced for 
the purpose of analyzing G2-compatible SH cosmologies. It is necessary to 
write out the Hubble-normalized evolution equations in component form (i.e. 
the 1+1+1+1 form!) rather than in the 1+1+2 form, so that we can impose a 
suitable spatial gauge condition. For convenience we introduce the following 
notation: 

iVc^ = 2iV+, N^s = ^[nI _^ ) > (4-24) 

and 

S12 = VSEg, Si3 = V3S2. (4.25) 

The Hubble-normalized equations in component form are obtained from the 
1+1+2 system ()2.47|) - ()2.63|) by expressing them in terms of the Hubble- 
normalized variables ()2.100|) - ()2.10H) . using ()2.84p - ()2.85|) . and are given be- 
low: 
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G2-compatible SH system (Hubble-normalized) 

Evolution equations for the gravitational field: 



EV = (g - 2)S+ + 3(S2^ + S3^) - 2{N'_ + N^) - ^^{vf - ^v') (4.26) 
E'_ = (g - 2)S_ + 2ANy, - 2i?Ex - 2N+N_ + ^/3{Es'^ - Ea^) 

+ f - «l) (4.27) 
E'^ = (g - 2)Ex - 2AAr_ + 2i?E_ - 2N+N^ + 2V3E3E2 





AT 7^ 
+ V 3— t'2t'3 


(4.28) 


S3' 


= (g - 2 - 3E+ - V3E_)E3 - {R + V3Ex)E2 + Vs^^viVi 


(4.29) 




= (g - 2 - 3E+ + V3E_)E2 + (i? - V3Ex)E3 + ^/3^v,V3 
= (g + 2E+)A^+ + 6E_A^_ + GE^A^x 


(4.30) 




(4.31) 




= (g + 2E+)Ar_ + 2Ar+E_ - 2RNy, 


(4.32) 


N' 


= (g + 2E+)iVx + 2A^+Ex + 2RN_ 


(4.33) 


A' 


= (g + 2E+)A, 


(4.34) 


where 








g = 2S2 + |G';i[(37 - 2)(1 - v^) + 27^;2]0 - fl^ 


(4.35) 




E^ = E^ + E^ + S^^ + Sa^ + E32 


(4.36) 




G± = 1 ± (7 - 1)^;2, . 


(4.37) 



Constraint equations: 

= (Cg) ^1- N1- Nl-A^ -E^ -Q-Qj^ (4.38) 

= (Cc)i = 2AE+ + 2(E_7Vx - SxiV-) + (4.39) 

= (Co) 2 = -(3A + V3Arx)E3 - {N+ - V3N^)^2 + V3]^V2 (4.40) 

= (Cc)3 = -(3^ - \/3iVx)E2 + (iV+ + \/3iV_)E3 + V3^^;3 • (4.41) 
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Evolution equations for the matter: 

= G-^[2G+q - (37 - 2) - (2 - 7)f2 - 7S„/3i;V + 2-fAvi]n (4.42) 

v[ = {M + 2S+)t;i - 2 V3(S3t;2 + ^2V3) + V^Ny, {vj - v^) + 2V3N.V2V3 - Av^ 

(4.43) 

v'^ = {M-^+- V3S_ + V3Ny,vi)v2 -[R+ VsSx - (iV+ - V3A^_)f i]f3 

(4.44) 



where 



(4.45) 
(4.46) 



M = G 



-1 



(37-4)(l-t;2) + (2-7)E,^t;V 
[G+ - 2(7 - l)]Av^ 



(4.47) 



= S+(t;2' + - 2vf) + V3J:4vI - vl) + 2V3Sxt;2t^3 

+ 2v^(S3i;2 + J:2V3)vi . (4.4J 



In order to derive the v'^ equations^ from the evolution equations ()2.62|) - 
()2.(iH|l for qa, we first obtain f2' and Q'^ by converting ()2.(il|) - ()2.(iH|) to Hubble- 
normahzed form. On differentiating Qa = 'jG'^'^^lva, we obtain 



(4.49) 



Contracting with and simphfying gives 



71; 



G 



2 7G_^ 



The resulting equation for f ^ is 

(^2)' = ^(1 - v') [(37 - - ^apv^v^ - 2(7 - lyAv,] 



(4.50) 



(4.51) 



^See van Elst & Uggla 1997 jTS,, page 2682. 
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which we use to evaluate the doG^ term in ()4.49|1 to obtain v'^. 

It follows from (ITT)T)|l - (pTT|l and that the Hubble-normalized Weyl 

curvature variables are given by 

£+ = 1(1 + + + Sa^) - |(E2 + El) 

+ m+N^_) + \^ivf-lv^) (4.52) 

£:_ = |(1-2S+)S_-^(S3^-S,^) 





9 1 AT , 9 AT AT , ,0 2\ 

1 AiVx + |iV+iV_ + 2^ (^2 ^3) 


(4.53) 




= |(1 — 2S_|_)Sx — -^S3S2 






+ |AiV_ + |iV+iVx + -^^t;2t;3 


(4.54) 


£3 


= 1(1 + S+ - v^S_)S3 - 75SXS2 - ^^^1^2 


(4.55) 


£2 


= 1(1 + E+ + V3S_)S2 - 75^x^3 - ^^^1^3 


(4.56) 




= -A^_S_ - A^xSx 


(4.57) 


n. 


= |ASx - |iV+S_ - S+Ar_ 


(4.58) 




= - |iV+Sx - S+iVx 


(4.59) 


^3 


= -|(A + V3iVx)S2 - ^iV_S3 + ^^^3 


(4.60) 




^l(A-v/3iVx)E3+;,iV_E2-^^^^2. 


(4.61) 



Note that the (Cc)2 and {Cc)^ constraints have been used to simplify 7-^2 and 
7-^3 respectively. 

Recall from Section that for SH cosmologies of Bianchi types VI and 
VII, the spatial curvature variables are related via 

A^ = h det Nab, det Nab = iV+ - 3(iV! + A^^^ ^ q_ ^^ q^) 

Thus the Bianchi VI and VII state space consists of "layers" of Bianchi VI/i 
and VII/i state spaces. 

For SH cosmologies with generic G2 action, a simple way to specify the 
spatial gauge is to use a time- dependent rotation of the spatial frame vec- 
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tors {e^} to make 62 parallel to a KVF. It then follows from Lemma 1 in 
Appendix [X] that 

N+ = V^N_, R = -VsSx . (4.63) 

We shall refer to this gauge as the Killing spatial gauge. This gauge, however, 
has the disadvantage that it cannot be used for Bianchi VII cosmologies, 
because it entails det A^^s < 0, as follows from (j4.62|) . 

A second possibility is to rotate the frame {e^} so that S3 = 0. The 
evolution equation ()4.29p and the constraint ()4.40|) then give 

Eg = 0, R= -VsSx + (A^+ - \^N_)vi . (4.64) 

We shall refer to this gauge as the shear spatial gauge. For the invariant set 
{v2 = 0, E2 7^ 0}, the (Cc)2 constraint dOOl) implies A^+ = V3N_, and 
simplifies to 

S3 = 0, A^+ = V3N_, R = -VsSx . (4.65) 

On the other hand, for the invariant set S2 = corresponding to SH cosmolo- 
gies with OT G2 action, the {€0)2 constraint ()4.4m) implies V2 = 0, leaving 
fOl)) the same. 

A third alternative is the tilt spatial gauge, where one rotates the frame 
{ga} so that f2 = 0. The V2 equation ()4.44j) then gives 

V2 = 0, R = -VsSx + (A^+ - VSN^)vi . (4.66) 

Observe that the tilt spatial gauge and the shear spatial gauge have the same 
expression for R. 

We shall use the shear spatial gauge ()4.64p rather than the tilt spatial 
gauge because it has the advantage of displaying the tilt degrees of freedom 
explicitly. With this choice, the Hubble-normalized variables are 

X= {^+,^.,^^,j:2,N+,N.,N^,A,n,vi,V2,V3,nA), (4.67) 

together with the parameter 7. There are 13 variables and 4 constraints, 
making the state space 9-dimensional. 



^The state space of Bianchi VI^, or VII/i cosmologies for each fixed h is 8-dimensional, 
taking into account the extra constraint H4.62|l . 
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Constraints 

The stability of the constraints is of great concern when doing numerical 
simulations. From Uggla et al. 2003 [721 equations (A2)-(A3)] we obtain the 



evolution equations for the constraints: 

(Cg)' = 2g(CG) - |A(Cc)i (4.68) 

(Cc); = (2g - 2 + 2S+)(Cc)i - 2V3S3(Cc)2 - 2v^S2(Cc)3 (4.69) 

(Cc)'2 = (2g - 2 - S+ - V3S_)(Cc)2 -{R+ v^Sx)(Cc)3 (4.70) 

(Cc)'3 = (2g-2-S+ + V3S_)(Cc)3 + (i?- v^Sx)(Cc)2 • (4.71) 



We see that the constraints are stable into the past near the Kasner circle 
(as the coefficients of the first terms are positive). 

4.3 Generalized Mixmaster dynamics 

We now predict under what circumstances Mixmaster dynamics can occur in 
(ji'2-compatible SH cosmologies. We use the gauge choice ()4.64|1 . but we shall 
show that the main result, Proposition 14.11 is independent of this choice. 

We begin by describing the Kasner equilibrium points. Within the state 
space of G2-compatible SH cosmologies, it is straightforward to verify that 
the Kasner circle /C, now given by 

(Sx,S2,iV+,iV_,iVx,A) =0, (4.72) 
+ = 1, (4.73) 

n = o = nA, va = o, (4.74) 

is a set of equilibrium points, and satisfies q = 2. In fact, by assuming 

n = = nA, {N+,N_,Ny,,A) =0, (4.75) 



one quickly arrives at = 
fjOH|l - (IOn|l imply 



1 and g = 2 through (p3H|l and ^i^B- Then 
S2 = = 0. (4.76) 
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where 



Evolution equations (jOH|) - (jOK|l and (p3Tjl with pTT^ and ^^^7^ give 

v[ = {M + 2i:+)vi (4.77) 

t;^ = (M - S+ - V3S_)t;2 (4.78) 

v'^ = (M - ^+ + V3^^)v3 (4.79) 

(^')' = c^(l - ^')i(37 - 4)t;^ - , (4.80) 

M = G:^[(37 - 4)(1 - f^) + (2 - 7)S„^f V] (4.81) 
S^^t;V = S+(t;^ + vl- 2vf) + V^^-ivj - vl) . (4.82) 

Equihbrium points with Va = give the Kasner circle JC. Equilibrium points 
with one non-zero tilt variable and f ^ = 1 give 

= (±1, 0, 0), (0, ±1, 0), (0, 0, ±1). (4.83) 

Hence we have 6 more Kasner circles )C±a of equilibrium points, given by 

(Sx,S2,iV+,iV_,iVx,A) = 0, (4.84) 
+ = 1, (4.85) 
n = = nA, = (±1,0,0), (0, ±1,0), (0,0, ±1). (4.86) 

The ]C±a equilibrium points are unphysical due to the extreme tilt f ^ = 1, 
but we shall see that they nevertheless play a significant role dynamically. 
Equilibrium points with one non-zero tilt variable and v"^ < 1 gives 

-i(37-4) if.;i^O 
v^S_ = 37-4 if W2 7^ (4.87) 
V3S_ = 37-4 if t;3 ^ 0. 

There are equilibrium points with two non-zero tilt variables, but they occur 
at Qa and and not on the entire Kasner circle. There are no equilibrium 
points with V1V2V3 7^ 0, since ()4.77|) - ()4.79|) require that 

M + 2S+ = = M - S+ - V3S_ = M - S+ ± VSS. , 
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which imphes S+ = S_ = 0, contradicting + = 1. The Kasner 
solutions are in fact characterized by ()4.75p . and thus have other (non- 
equihbrium) representations, which we shall encounter very soon in this chap- 
ter. 

In order to proceed with the analysis we first restrict our consideration 
to solutions that satisfy 

lim = 0, lim = 0, (4.88) 

t—*—oo t— » — oo 

based on the BKL conjecture I (see Section lOj) . Second, based on the Kasner 
Attractivity Conjecture below, we restrict our considerations to solutions 
whose orbit enter an arbitrarily small neighbourhood of one of the Kasner 
circles. Numerical simulations suggest that this property is satisfied, and 
hence we formalize it ClS db conjecture: 

Conjecture 4.1 (Kasner Attractivity Conjecture). The orbits of all 
G2-compatible SH cosmologies, except for a set of measure zero, enter an 
arbitrarily small neighbourhood of one of the Kasner circles. 

Consider linearizing the evolution equations ()4.28|1 . ()4.30|) - ()4.34|) for the 
spatial curvature and off-diagonal shear variables ()4.72|1 on the above Kasner 
circles: ^ 





= -2V3S_Sx 


(4.89) 




= (-3S+ + y3S_)S2 


(4.90) 




= 2(1 + S+ + V3S_)iV22 


(4.91) 




= 2(1 + S+ - V3S_)iV33 


(4.92) 


N' 


= 2(1 + j:+)n^ 


(4.93) 


A 


= 2(1 + S+)A, 


(4.94) 



where (S_|_, E_) = (cos0, sin0), and is a constant. It follows that the first 
four variables are unstable on some arc of the Kasner circles, i.e. will act as 
triggers, while A^x and A are stable on the Kasner circles, into the past. It 
is plausible that the stable variables tend to zero: 

lim (A^x,^) = 0, (4.95) 
t— >— 00 

^Recall from that N22 = N+ + V^N^, N33 = iV+ - V3N^. 
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and this conclusion is strongly supported by numerical simulations. 

We recognize that N22 and N^s are two of the curvature triggers intro- 
duced in Section l^m In the present situation, however, the curvature transi- 
tion sets 7/V22 and T/Vga also describe transitions on the Kasner circles /C±q,, in 
addition to /C. The new triggers here are Ex and S2: the off-diagonal shear 
components. The Ex transition orbits lie on the sphere -|- -|- = 1. 
Evaluating ()4.26|) on the sphere gives S'^ = 0, which implies 



on each Ex transition orbit. Similarly, the S2 transition orbits lie on the 
sphere + = 1. Evaluating K^ - KTfl on the sphere gives 

(S+ -|- -\/3S_)' = 0, which implies 



on each S2 transition orbit. The Ex and S2 transition sets are shown in 
Figure By comparing Figure with the dashed lines in Figure ITTk . we 
see that these transition sets connect different representations of the same 
Kasner solutions. These transitions are in fact an artifact of a rotating spatial 
frame (non-zero fi^), and the transition sets actually represent the Kasner 
solutions. Therefore we shall call the off-diagonal shear components the frame 
triggers, with corresponding frame transition sets, denoted and • ^ 

The triggers that lead to the instability of the various Kasner arcs are 
shown in Figure 14. 4| which is obtained using Figures I4.2H4.3I The essential 
point is that 

Mixmaster dynamics occurs if and only if each Kasner arc has 
at least one trigger, and so if an arc has no triggers, then it 
will form part of the past attractor, and Mixmaster dynamics 
will not occur. 

Table mU lists the variables for the subclasses with generic G2 action and 
the additional restrictions. 



^This artifact is beneficial, ff one chooses a Fermi-propagated frame, i.e. fla = 0, then 
the Kasner equihbrium points form a hypersphere 



S+ = const. 



(4.96) 



= const. 



(4.97) 



= 1 



instead of a circle + E'i 



= 1. 
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Figure 4.3: The frame transition sets corresponding to the triggers S 
S2, respectively. The arrows indicate evolution into the past. 



Sx,iV33 




^2? N22 



Figure 4.4: The unstable Kasner arcs and their triggers. 



CHAPTER 4. SH DYNAMICS 



56 



Table 4.1: The Hubble- normalized variables for the G'2-compatible SH cos- 
mologies with A = and generic G2 action (see Figure I3.3|) . relative to 
the shear spatial gauge. Trigger variables appear in bold. Recall that 



Subclass State vector X Additional restriction 



With 3 tilt de^ 


rrees 


of freedom: 








VI, VII (S4 




,Ey,,i:2,N+,N_,Ny,,A,n,vi,V2,V3) 


none 






IV (S^ 




,Ex,'^2,N+,N^, Ny.,A,n, vi,V2,V3) 


_ 


3{N^ + A^2 


) = 


With 2 tilt dej 


rrees 


of freedom: ^2 = 


N+ = 






VI (S^ 




,-E^,-E2,N_,N^,A,n,vi,V3) 


A^x = 


-^/3A if /i 


1 

~ 9 


IV (S4 




,'Sx,'^2,N-,A,n,vi,V3) 


iVx = 







V (S4 




,T;x,T,2,A,n,vi,V3) 


A^+ = 


A^_ = A^x = 


= 


With 1 tilt dej 


rree of freedom: t;2 = f 3 = 


A^+ = 


V^N_, A^x 


= V3A 


vr_i (s^ 






none 






II ' (S^ 






A^x = 


^ = 




Non-tilted: 


= 




Ar+ = 


V3Af-, A^x 


= V3A 


VP 1 (S^ 






none 







Figure 14.41 and Table 14.11 provide the basis for the following proposition. 

Proposition 4.1. Consider the class of perfect fluid G 2- compatible SH cos- 
mologies which satisfy the Kasner Attractivity Conjecture. 

i) If the G2 action is generic, then Mixmaster dynamics occur if and only 
if the anisotropic spatial curvature Sap is non-zero. 

a) If the G2 action is non-generic, then Mixmaster dynamics does not 
occur. 

Proof. From Figure 14.41 and Table 14.11 we see that for models of all Bianchi 
types with a generic G2 action, each arc on the Kasner circles has at least 
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one trigger, unless = = A^x = (Bianchi type V, where Sap = 0). It 
follows that Mixmaster dynamics occur if and only if Sajs 7^ 0. 

On the other hand, if the G2 action is non-generic, the state vector X in 
fl4.67|) is restricted as follows. For SH cosmologies with OT G2 action, from 
(|3.5|) we have 

S2 = 0. (4.98) 
For SH cosmologies with one HO KVF, from ()3.2|) we have 

Ar+ = = = Sx. (4.99) 

For SH cosmologies with diagonal G2 action, we have 

S2 = 0, N+ = =0 = ^^. (4.100) 

Figure 113] shows that in each case there is an arc with no triggers, implying 
that Mixmaster dynamics does not occur. 

□ 

Comment: To show that this result is independent of using a particu- 
lar gauge, namely the shear spatial gauge ()4.64|) . we note that the infinite 
sequences of transitions that arise lead to infinite sequences of changes in 
physically significant scalars, for example, the scalars formed from the Weyl 
curvature tensor, as we describe later. Since these changes are independent 
of the choice of gauge, the occurrence of Mixmaster dynamics is independent 
of the choice of gauge. However, the choice of spatial gauge does affect which 
triggers occur and hence the details of the dynamics.^ For example, in the 
Killing spatial gauge ()4.63|1 . the triggers are Sx , S2 and S3. It should be 
kept in mind that the frame transitions are merely an artifact of the rotating 
spatial frame, and do not lead to a change of Kasner state. ^ 



^The real, physical dynamics, and artifacts of the rotating spatial frame both affect 
the triggers and other gravitational field variables. To separate these effects, we observe 
rotation-invariant scalars formed these variables (see (|4.117|l - (|4.119|l below). 

®As an extreme example, an artificial type of Mixmaster dynamics consisting of an 
endless sequence of frame transitions occurs for models of Bianchi type V if the Killing 
spatial gauge (|4.63|) is used, since the only triggers are Sx, S2 and S3. The rotational 
freedom can be used to set E3 = to eliminate this artificial Mixmaster dynamics (the 
arc {T1Q2) becomes stable). 
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In addition to ()4.95j] . there appear to be other hmits that help to charac- 
terize the asymptotic dynamics. From the hnearized equations ()4.89p - ()4.92p 
we obtain 

(iV22iV33)' = 4(1 + S+)(iV22iV33) (4.101) 
(iV22Sx)' = 2(1 + S+)(iV22Sx) (4.102) 
(iV33S2)' = (2 - S+ - V^S_)(iV33S2) . (4.103) 

It follows that these products are stable on the Kasner circles into the past, 
suggesting that they tend to zero as t — — oo: 

lim (iV22iV33,iV22Sx,iV33S2) = . (4.104) 
t— »— oo 

This claim is strongly supported by numerical simulations. These limits can 
be interpreted as stating that the two triggers in each product cannot be 
active simultaneously. 

Multiple transitions 

Returning to Figure 14. 4| notice that there is more than one trigger on some 
Kasner arcs, which means that two triggers can become active simultane- 
ously, leading to more complicated transitions between Kasner states. These 
so-called multiple transitions were first studied in detail by Hewitt et al. 
2003 ^36j in the context of non-tilted exceptional Bianchi Vl^^^g cosmologies. 
The orbit of a multiple transition is not a straight line when projected on the 
(S+, E_) plane. See Figures 3 and 4 in [20] • Two of the three possible multi- 
ple transitions involve the pairs (A''22, S2) and (A''33, Ex), which describe the 
Taub vacuum Bianchi II solutions relative to a non-Fermi-propagated frame. 
The third one, involving the pair (E2,Sx), describes the Kasner solutions 
relative to a non-Fermi-propagated frame. The possibility of multiple transi- 
tions is a significant difference from the "standard" Mixmaster dynamics (of 
Bianchi VIII and IX). On the basis of detailed numerical simulations, Hewitt 
et al. conjecture that the probability of multiple transitions occurring tends 
to zero as t — ^ —00. The numerical experiments suggest that multiple tran- 
sitions keep occurring as t ^ — 00, but that their strength, as measured by 
|A^22E2|, |A'^33Sx| and |S2Sx|, keeps decreasing, i.e. that 



lim (iV22S2,iV33Sx,S2Sx) =0 



(4.105) 
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We caution that this result is not very conclusive, and we shall revisit multiple 
transitions in the context of G2 dynamics in Chapters [7| and El 

Tilt 

What happens to the tilt variables Wq,? To determine the behaviour of Va, 
linearize on IC, and also and (ICTll on IC±i: 

on/C: v[ = - A + 2J:+)vi (4.106) 

(37 - 4 - S+ - V3^-)v2 (4. 107) 

(37-4-S+ + V3S_)t;3 (4.108) 



^2 



on /C±i: v'2 = (-3S+ - V3S_)t;2 (4.109) 

fg = (-3S+ + V3J:^)v3 (4.110) 

(1 - v'y = -2^— +(1 - v') . (4.111) 

2-7 

For the purpose of this thesis, we shall focus on models with only one tilt 
degree of freedom, i.e. Vi ^ 0, V2 = V3 = 0. Equations ()4.106|) and ()4.111|] 
give the linearized solutions 

V, = (t;i)oe(=^^-^+2^+)* (4.112) 

(37-4+2E+) 



Vl-v^ = {Vl-v^)oe . (4.113) 

Hence, the tilt is unstable into the past on the following arcs: 

/C arc(E+ < -|(37-4)), /C±i arc(S+ > -|(37 - 4)). (4.114) 

This creates transition orbits between /C and /C±i, as illustrated in Figure 117^ 
The unstable arcs are highlighted, and the arrows indicate evolution into the 
past. For convenience we augment Figure H^ bv adding the unstable arcs for 
f 1 to produce Figure 14.61 There are analogous transition orbits between /C 
and /C±2, and between K. and /C±3. There are also transition orbits between 
the extreme tilt Kasner circles. ^ 



^We refer the reader to the work by Uggla, van Elst, Wainwright & EUis 2003 in 
the context of Gq cosmologies. 
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Figure 4.5: The tilt transition for Vi. 



Comments: The nature of the tih transition means we cannot make non- 
Mixmaster dynamics into Mixmaster by adding only tilt transitions. However 
adding tilt degrees of freedom to a non-tilted model also adds off-diagonal 
shear variables and it is these variables that lead to the occurrence of Mixmas- 
ter dynamics. From Table l^?Tl we see that the tilted Bianchi II cosmologies 
form the prototype for Mixmaster dynamics with tilt. We refer to the anal- 
ysis of tilted Bianchi II cosmologies by Hewitt et al. 2001 The class 



^"^The spatial gauge used in this paper is a permutation of the one used in this thesis. 
The indices they use are "one lower" than ours, e.g. they have 

iVn^O, V3^0, S3i=0, (4.115) 

while we have 

N22 7^ 0, VI ^ 0, S12 = 0. (4.116) 
They define T,± relative to their 23-subspace, as we do in H4.24|) . 
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Figure 4.6: The unstable arcs on IC and IC±i and their triggers. 

of tilted exceptional Bianchi Vl^j^^g cosmologies has not been analyzed, but 
we can now predict, without additional work, that its Mixmaster dynamics 
is that of the tilted Bianchi II cosmologies. 

Physical aspects of Mixmaster dynamics 

The physical aspects of Mixmaster dynamics are described by rotation- invariant 
scalars, such as the shear scalar E and the electric and magnetic Weyl scalars 
S and H: 

E^^E^ + E^ + E^^ +E2^ + E32 (4.117) 
S^^£l + £^_ + Sl+ E^"" + (4.118) 

= nl + nl + nl+n3^ + n2'^. (4.119) 

Furthermore, the scalar S'^ — Ti^ is boost-invariant. The shear scalar S is 
close to 1 except during curvature transitions. The other invariant of E^/j is 
the determinant det E^^, which is also equal to ^Sa/^E'^'^'E^". det Eq,^ is close 
to constant except during curvature transitions. The electric Weyl scalar £ 
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is equal to a/2 — det when = 1. It is close to constant except during 
curvature transitions. The magnetic Weyl scalar Ti is close to zero except 
during curvature transitions. We refer to WE, Section 11.4 for the graphs 
of these scalars under the standard Mixmaster dynamics. There are other 
rotation-invariant scalars, such as Q^, SapS"^ and detiS^^. The principal 
feature of Mixmaster dynamics is that the Hubble-normalized Weyl scalars 
are bounded but do not have a limit in the past asymptotic regime. 

Another physical aspect is the fluid vorticity. Although the timelike con- 
gruence eo is vorticity-free, the fluid congruence is not necessarily so. From 
Appendix C of Hewitt et al. 2001 ^25,, the fluid vorticity for G2-compatible 
SH cosmologies is given by 

^1 = W^T^2^''''''cVdVi (4.120) 
Zn 1 — f ^ 



Wa = — 
2B 



{Na^ + eA^A)vc + -N^^vcVdVa 

1 — 



(4.121) 



where 

1 - \{v^ + T.^pv'^v^ + 2Avi] 



B 



2 



Observe that the fluid vorticity is zero if f2 = = 0. 

For the case V2 = 0, fs 7^ 0, the spatial gauge is given by ()4.65|) . which 
implies N'^^vc'Vd = 0. It follows that the only non-zero vorticity component 
is 

W2 = ^iV3N^ + A)v3 . (4.122) 
ZJd 

The limit ()4.95|) then implies that W2 tends to zero. 

For the case with three tilt degrees of freedom, since the limit of N^^vcVd 
does not exist, the limit of the fluid vorticity does not exist a.s t —00, i.e. 
Mixmaster dynamics is reflected in the fluid vorticity. 

Thus we have the following proposition. 

Proposition 4.2. Consider G 2- compatible SH cosmologies with generic G2 
action. 

• If there is one tilt degree of freedom (i.e. va = 0), then the fluid 
vorticity is zero. 

• If there are two tilt degrees of freedom, then the fluid vorticity tends to 
zero as t —>■ —00. 
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• // there are three tilt degrees of freedom, then the limit of the fluid 
vorticity does not exist as t — oo. 

We revisit Figure and include the total number of curvature and frame 
triggers, the dimension of the state space, and the asymptotic behaviour 
of the fluid vorticity (see Figure 14. 7|) . Note that for Bianchi VI and VII 
cosmologies, the dimension listed is for the state space with a fixed parameter 
h. For A > 0, increase the dimension by one. 

Before giving a summary we note that there is some uncertainty as regards 
the behaviour of the tilt as t — > — oo. The analysis of the triggers and of the 
tilt transitions suggests that along a typical infinite sequence of transitions, 
the tilt scalar v will not approach a limit: 

lim V does not exist. (4.123) 

t—^—oo 

However, while performing numerical simulations of Bianchi II models with 
7 = |, Bridson 1997 [TS] discovered that, despite the existence of tilt tran- 
sitions from K,±i to /C that transform v from ±1 to 0, v tends to get stuck 
at ±1. He argued that, due to the arrangement of the transition orbits, the 
state vector X spends more time close to the part of lC±i where 1 — f ^ is 
stable and hence 1 — f ^ is forced even closer to zero, so that when it has an 
opportunity to activate, i.e. when X is close to the unstable arc on lC±i for 
1 — f ^, it needs more time to do so and loses out to the curvature and frame 
triggers. This tendency will, however, depend significantly on the value of 7. 
If 7 is close to 2, then v should have a strong tendency to stay close to zero, 
for two reasons: (1 — v"^) is unstable on most of /C-ti, and also has a high 
growth rate into the past due to the (2 — 7) factor in the denominator (see 
equation (j4.113|) ). Hence one expects that somewhere between 7 = | and 
7 = 2 there is a critical value of 7 where one tendency replaces the other. 
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/ciet uab 7^ 0, 


= hdet 


riAB 




det uab 


= 




VI, VII 

7D 


4 triggers 






IV 


4 triggers 




Stilt 


> 


7D 


Stilt 



lim W does not exist 

t^—oo 



uab = 





3 triggers 


VI 




6D 


2 tilt 



IV 



6D 



3 triggers 
2 tilt 





2 triggers 


V 


5D 


2 tilt 



VI* 1 



6D 



VP 1 



5D 



W^^O, lim iy = o 

t— >— oo 



3 ttiggers 
Itilt 





3 triggers 


II 




5D 


Itilt 



3 triggers 

tilt 



W 



Yj* 3 triggers 

4D 9 vacuum 



Figure 4.7: SH cosmologies with generic G2 action, indicating the Bianchi 
group type, the dimension of the Hubble-normalized state space, the number 
of triggers and the number of tilt degrees of freedom. Here W is the Hubble- 
normalized fluid vorticity scalar, defined by W"^ — + WaW^. 
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Summary 

To recapitulate, the generalized Mixmaster dynamics that we have discussed 
in this section can be described by the following limits: 

BKL conjecture I: lim = 0, hm VLp, = 0, (4.124) 

Stable variables: lim (A, A^x) = 0, (4.125) 

t— >— oo 

Non-overlapping triggers: lim (A^22A^33, ^^22^x5 A^33S2) = 0, (4.126) 

t— >— 00 

Trigger variables: lim (Ai'22, -^^33, S2, Sx) do not exist. (4.127) 

t^— 00 

Unlike the case of standard Mixmaster dynamics discussed in Section li?H 
to date no proofs have been given concerning generalized Mixmaster dynam- 
ics (e.g. the limits ()4.124j) - ()4.126|l which serve to restrict the past attractor). 
The non-tilted Bianchi Vl^^^g solutions and the tilted Bianchi II solutions are 
prototypes for the new features of generalized Mixmaster dynamics. Under- 
standing these two special cases, and having our formulation of the evolution 
equations for the whole hierarchy, gives qualitative information about the dy- 
namics of the whole class and provides equations for doing detailed analyses, 
and perhaps for providing proofs of some of the above limits. 



4.4 Non-Mixmaster dynamics 

In this section, we discuss the past asymptotic dynamics of SH subclasses 
with non-generic G2 action. The conditions for OT G2 action, HO KVF and 
diagonal G2 action are given by ()4.98p - ()4.100|) respectively. Figures I3.4H3.5I 
summarized the hierarchy for these SH subclasses. Tables I4.2H4.4I list the 
variables and the additional restrictions. 

For each entry in Tables I4.2H4.41 we combine the unstable arcs of the 
curvature, frame and tilt transition sets in Figure 14.61 and look for the re- 
maining stable arcs on the Kasner circles K. and /C±i. In each case, sequences 
of transition sets will terminate at some stable arc. For each class the past 
attractor is an arc or union of arcs on the Kasner circles. The result is that 
each orbit, except for a set of measure zero, is past asymptotic to a Kasner 
equilibrium point. 

For SH cosmologies with OT G2 action ()4.98|) . the S2 trigger is zero. 
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Table 4.2: The Hubble-normalized variables for the G2-compatible SH cos- 
mologies with A = and OT G2 action (see Figure . Triggers appear in 
bold. 



Subclass State vector X Additional restriction 



With 1 tilt deg 


;ree 


of freedom: ^2 = ^3 = 


S2 = 











VII (S+ 




-,'£y,,N+,N_,Ny,,A,n,vi) 


none 










VI (S+ 




^,'£^,N.,N^,A,n,vi) 


N+ = 










IV (S+ 


,s_ 


.,-£y,,N_,A,Q,vi) 


N+ = 


■ VSN^, 




= 




II (S+ 






N+ = 


■ V^N^, 




= = 


A 


Non-tilted: Va 


= 




S2 = 











VII (S+ 




-,'S^,N+,N-,N^,A,n) 


none 










VI (S+ 




^,j:^,N.,N^,A,n) 


N+ = 










IV (S+ 




,^^,N.,A,n) 


N+ = 


■ VSN^, 




= 




II (S+ 






N+ = 




A^x 


= = 





making the Kasner arc (T3Q1) on /C in Figure stable. The corresponding 
arc on ]C±i is unstable due to the tilt. Thus the conjectured past attractor 
is given by 

A' = IC arc(T3Qi). (4.128) 

The past attractor for the non-tilted case (i.e. Vi = 0) is also given by ()4.128|) . 
except for Bianchi II, whose past attractor is a bigger set 

A- = IC ^iciT^Q^T^QiT^). (4.129) 

For SH cosmologies with one HO KVF (|4.99p . there is only one trigger, S2. 
With one tilt degree of freedom, from Figure |33l we see that the conjectured 
past attractor is given by 

A- = [/C arc(S+ > -|(37 - 4)) n /C arc(Q2TiQ3T2) ] 

U [/C±i arc(S+ < -i(37-4)) n /C±i arc(g2Tig3T2) ] . (4.130) 

For non-tilted models, the past attractor is the point (S+, S_) = (— |) on 
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Table 4.3: The Hubble-normalized variables for the (j'2-compatible SH cos- 
mologies with A = and one HO KVF (see Figure 13. 5|) . S2 is the only 
trigger. 



Subclass State vector X 


Additional restriction 






With 2 tilt degrees of freedom: ^2 = 
VI (S+,S_,E2,iVx,A^^,^^i,^'3) 

V (i;+,s„,S2,^,r2,fi,f3) 


N+ = VSN- = = Sx 
Ny, = -VSA if /i = -5 
A^x =0 






With 1 tilt degree of freedom: ^2 = f 3 = 
VI% (S+,S_,S2,Af^,^^i) 


A^+ = VSN^ = = Sx, 
none 


A^x 


= V3A 


Non-tilted: Va = 

VI* 1 (s„,S2,Af^) 


A^+ = V3N_ = = Sx, 
none 


A^x 


= V3A, 



KL. Finally, we note that the past attractor for the case with two tilt degrees 
of freedom is also given by ()4.130|) . 

For SH cosmologies with diagonal G2 action (j4.1(Jm) . there are no triggers. 
From Figure 1^7^ we see that for tilted models (fi 7^ 0) the conjectured past 
attractor is given by 

A~ = K, arc(S+ > -i(37 - 4)) U /C±i arc(S+ < -i(37 - 4)). (4.131) 

-'^-'^It turns out that the set for stable arcs on /C±3, given by 

/C±3 arc(S]+ - V3I]_ > 87 - 4) n /C±3 aTc{QzT2QiT:i) n /C±3 arc(Q2riQ3T2) , 

is empty, and the arcs for which U3 = is stable on lC±i is identical to K,±i arc(Q2T'iQ372) 
(compare (|4.1U8|) with H4.90|l ). Finally, the stable arcs on /C, given by 

K. arc(S+ > -^(87 - 4)) n /C arc(I]+ - V3I]_ < 87 - 4) n /C arc(g2TiQ3r2) , 

simplifies to 

/C arc(S+ > -1(87-4)) n /C arc(g2TiQ3T2) . 
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Table 4.4: The Hubble-normalized variables for the G2-compatible SH cos- 
mologies with A = and diagonal G2 action (see Figure IH3j) . 



Subclass State vector X Additional restriction 



With 1 tilt de^ 


^ree of freedom: V2 = V3 = 


S2 = 


0, N+ = 


^iV„ = = Sx 


VI 


(S4 




none 






V 


(S4 




A^x = 


-- 




V (LRS) 


(S4 


-,A,Q,vi) 


A^x = 


: = 




Non-tilted: 


Va 


= 


S2 = 


0, N+ = 


y3A^_ = = Sx 


VI 


(S4 




none 






V 


(s_ 


-,A,n) 


A^x = 


: = S + 




I 


(5^4 




A^x = 


: = A 




I (LRS) 


(S4 




A^x = 


--0 = A = 





With zero tilt, the past attractor is given by 

A- = IC , (4.132) 
except for Bianchi V, whose past attractor is a smaller set 

A- = {S+ = 0, S„ = ±1}. (4.133) 

4.5 Future asymptotics 

In this section we discuss the future asymptotic dynamics of G'2-compatible 
SH cosmologies. 

The role of A 

The presence of a cosmological constant A > simplifies the future asymp- 
totic dynamics of SH cosmologies dramatically, as shown by Wald 1983 jHU] . 
He proved, for all SH cosmologies excluding Bianchi type IX, subject to cer- 
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tain restrictions on the matter content, that 

hm H = a/4, hm = 0, (4.134) 

t— >oo V t— >oo 

where H is the Hubble scalar and aa/3 is the shear tensor. In other words, 
the cosmological models are asymptotic to the de Sitter solution. 

Starobinskii 1983 [HZI discussed the asymptotic dynamics of cosmologies 
without symmetry assumptions and provided a conjecture for behaviour of 
the tilt variables, an issue that Wald had not addressed. He called the positive 
cosmological constant "the best isotropizer" . 

Our goal is to describe the effect of including a positive cosmological 
constant A from the perspective of the Hubble-normalized state space. First, 
including A increases the dimension of the Hubble-normalized state space 
by one, by adjoining the Hubble-normalized variable Q^. Second, the Gauss 
constraint (|4.38|) implies 

< < 1. (4.135) 

The condition 

= 1 (4.136) 

with and implies that 

n = 0, = 0, = = A^x = ^, q = -l, (4.137) 

and the evolution equation (j4.46|) for I^a is identically satisfied. Equations 
(j4.136|) - (j4.137|) thus define an invariant set. The evolution equation (j4.31|) 
evaluated on the invariant set ()4.136|) - ()4.137|) gives 

A^; = -A^+ , (4.138) 

which implies that equilibrium points must satisfy 

A^+ = 0. (4.139) 



Evolution equations ()4.43|) - ()4.45|) and ()4.51|) evaluated on this invariant set 
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give 



< = ^^^(l-^>a (4.140) 



{v'y = ^^^^il-v')v\ (4.141) 



which imply that the equihbrium points are 



Va = or v'^ = 1 if 7 7^ I 

Va = const., < f ^ < 1 if 7 = |. 



(4.142) 



We shall refer to the set (im?)|) - (l4.i;-i7jl . (I4.i;-i9;i and as the de Sitter 

invariant set. 

The evolution equation (j4.46|) for ^a, (I4.35p and (j4.135p imply that 

fi'A = (g + 1)^^A > , (4.143) 
which implies that Q\ is monotone increasing, and 

lim = 1, hm ^a = 0. (4.144) 

The presence of ^a thus has no effect on the past asymptotic dynamics, but 
destabilizes all equilibrium points on the set Q\ = into the future. That 
lim ^A = 1 implies 



t^oo 



limfi = 0, lim = 0, lim (A^+, A^_, A^x, ^) = 0. (4.145) 

As a result, the Weyl curvature variables also tend to zero: 

lim Sap = 0, lim Hap = 0. (4.146) 

Since there are no equilibrium points with < Q\ < 1, as t ^ oo the 
solutions must approach the de Sitter invariant set. The tilt variables have 
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the following limits: 




= if 7 < I 

= const. if 7 = I 

= const., lim f ^ = 1 if | < 7 < 2. 



t— >oo 



(4.147) 



See Lim et al. 2004 |1H], Appendix B for a proof. 
The case A = 

The future dynamics in the case A = is much more complicated than in the 
case A > 0, and depends significantly on the Bianchi type and on the equation 
of state parameter 7. The future dynamics for the non-tilted case has been 
completely analyzed, but much remains to be done as regards the tilted 
case. Based on the known cases - the Bianchi II and VIq subclasses (see jHS] 
and !32| respectively), non-tilted equilibrium points in the interior {VL > 0) 
of the state space, which act as the future attractor for non-tilted models, 
are destabilized by the tilt as 7 increases. The resulting bifurcations lead to 
the creation of tilted equilibrium points that act as the future attractor. 

Table 14.51 gives an overview of the work on tilted G2-compatible SH cos- 
mologies (all with A = 0) using Hubble-normalized variables. Hervik and 
Coley 2004 are currently analyzing the future dynamics of tilted Bianchi 
Vlh cosmologies, and a prehminary result suggests that a complicated, peri- 
odic orbit attractor exists for a certain range of 7. 

The class of SH cosmologies of Bianchi type VII is of particular interest 
because it contains the flat and open FL cosmologies as isotropic special cases 
(fiat FL as a special case of VHq and open FL as a special case of Vllh)- 

The Bianchi VHq state space is unbounded: is not bounded by the 
Gauss constraint ()4.38p . This is not an artifact of a bad choice of spatial 
frame, since the Weyl curvature scalars become unbounded as Nj^ 00 (A^+ 
appears in and p3H|l - (ji3n|) ). In other words, the Hubble- 

normalized Weyl curvature scalars are unbounded on the Bianchi VHq state 
space. The detailed analysis of the non-tilted Bianchi VHq cosmologies shows 



i^See WE, Wainwright, Hancock & Uggla 1999 "TS", Nilsson, Hancock & Wainwright 
2000 Horwood et al. 2003 Mi and Hewitt et al. 2003 36 . 

^■^See WE for references to work on tilted SH cosmologies using other approaches. 
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Table 4.5: Previous work on G2-compatible SH cosmologies with tilt and 
A = 0. 



Bianchi type 


G2 action 


paper 


Bianchi V 


generic G2 


Harnett 1996 jHO] 




HO KVF 


Harnett 1996 jHI] 




tilted diagonal 


Hewitt & Wainwright 1992 jSZj 




LRS tilted diagonal 


Collins & Elhs 1979 19j 


Bianchi II 


generic G2 


Hewitt et al. 2001 ^ 


Bianchi VIq 


generic G2 


Hervik 2004 p] 



that if 7 < |, they isotropize as regards the shear (i.e. T,± as r 00) 
even though A''+ 00, but that the Weyl curvature diverges if 7 > 1 [ZH! 
|56j . Preliminary investigation of the tilted models |17] indicates that no 
tilted equilibrium points exist for Bianchi VHq cosmologies. We expect that 
if 7 < |, they isotropize as regards the shear. This matter requires further 
investigation. 



Chapter 5 



The state space for G2 
cosmologies 

In this chapter, we use the 1+1+2 orthonormal frame formahsm as developed 
in Section IT^ to formulate the evolution equations for G2 cosmologies. We 
also introduce the notion of asymptotic silence and the silent boundary, which 
provides the link between G2 dynamics and SH dynamics. 

5.1 The evolution equations 

The starting point for deriving the evolution equations for G2 cosmologies 
is the 1+1+2 form of the orthonormal frame equations. As discussed in 
Section 12.31 we normalize the variables for G2 cosmologies using the area 
expansion rate /3, as motivated by van Elst et al. 2002 [73]. /^-normalization 
is advantageous only when used in conjunction with the so-called separable 
area gauge, introduced by van Elst et al. 2002 jTlj in the context of OT G2 
cosmologies. 

The separable area gauge condition is given by 

U = r. (5.1) 

By using the {Cu)i constraint ()2.77|) we find that A/", introduced in ()2.107|) . 
is independent of x, which implies that we may reparametrize t to set 

J\f=l. (5.2) 
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A special case of the separable area gauge is the timelike area gauge, for 
which 

U = r, Af=l, A = 0. (5.3) 

This gauge is equivalent to the one that has been widely used in the study 
of vacuum OT G2 cosmologies using the metric approach (see, for example, 
Rendall & Weaver 2001 [02] )• Although the timelike area gauge is simpler 
than the separable area gauge in a number of respects (see Sections 18. ![ 18.21 
and l8.3p . we shall resist setting A = 0, since the separable area gauge enables 
one to include all G2-compatible SH cosmologies as special will be 

explained in Section 15121 

As in Section 1^21 for convenience we introduce the following notation: 

Nc"^ = Nab = ^{^n^ _^ ) > (5-4) 

and 

E12 = \^S3, Si3 = V3S2. (5.5) 

The 1+1+1+1 /^-normalized evolution equations and constraints then follow 
from the system ()2.47|l - ()2.(i3|l . and are given below. In the process, we replace 
H hj (3 using the relation 

H={l-E+)f3 (5.6) 

(see fl2.103|) ). and normalize the variables using ()2.100|) - ()2.101|) . The vari- 
ables q and r are introduced and are given by ()2.114p and ()2.115|) respectively. 
We remind the reader that 



(5.7) 
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G2 system (/^-normalized) 

Evolution equations for the gravitational field: 

dtEi^ = {q + S^+)Ei^ (5.8) 

dtA = {q + 3S+)A (5.9) 
dt^+ = -{q + 3S+)(1 - S+) + 2(S+ + S?. + S2 + + Sg^) 
+ ^G^^liSj - 2) + (2 - jy]n - - |(di - 2A)r 

d^N^ = {q + 3E+)7V+ + 6E_A^_ + 6Ex A^x - diR 
dtT,- + diNx = {q + 3E+ - 2)E_ + 2AiVx - 2EEx - 

+ V3(E3^-S2^) + ^£(.i-.i) 

dtN^ + di^ = (g + 3S+)iVx + 2iV+Sx + 2RN^ 
dtT.^ - diN_ = {q + 3S+ - 2)Sx - 2AN_ + 2i?S_ - 2N+N^ 

G+ 

dtN^ - 5iEx = (g + 3E+)A^_ + 2A^+E_ - 2RNy, 

dtEs = {q-2- V3E_)E3 - (i? + V3Ex )E2 + Vs^^viVi 

G+ 

dt^2 = {q-2 + V3E_)E2 + iR- V3Ex)E3 + V3^t;ii;3 , 



where 



1 + 3(-E2^ - Sa^ + E^ + E^ + Ar2 + Ar2) _ ^2 _^ 2Ar 
+ 3G+i[(7 - 1)(1 - i;2) + jvfjn - SQa 
r = -3^E+ - 3(iVxE_ - A^.Sx) - 

2 C;r_|_ 

G± = i±{^-iy, v^ = vf+vi+vi. 

Constraint equations: 

= (Cg) = 1 - 2E+ - E2 - E2 - E2^ -^s^-Q-Qa 
-NI-NI + l{di - r)A - A^ 

= (Cc)2 = {di - r)E3 - {3A + V3iVx)E3 - {N+ - \/3A^_)E2 +V3]^V2 

G+ 

=(Cc)3 = (^1 - r)E2 - (3^ - V3iVx)E2 + (iV+ + ^7V_)E3 +Vs^V3 
= {Cp) = {di - 2r)^K . 
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Evolution equations for the matter: 



dt^ = -—diQ + — 

- - + Vidi lnG+ + 2v4i;i 



2^(g+l)-3(l-S+)(l + lt;^^ 



(5.25) 



Sft'i = —vidiVi + ^1 lnG_| 



7-1.. ^,2 



7 



+ (M + 2S+)i;i - 2y3(S3i;2 + ^2V3) 

V2V3 — Av 

dtV2 = -vidiV2 + (M - S+ - V3S_ + V?,N^vi)v2 
-[R+ VsSx - (iV+ - V?>N^)vi]v3 

dtV3 = -v^divs + (M - S+ + V3S_ - V3Ny,vi)vs 

+ [R- v^Sx - (iV+ + V3Ar_)t;i]t;2 
St^A = 2(g + l)fiA , 



(5.26) 

(5.27) 

(5.28) 
(5.29) 



where 



M = G-^ 



+ 



(7 - 1)(1 - v^)divi - (2 - -r)vidi In Gh 



7 



-(2 - 7)(1 - v^)vi{di InQ- 2r) + G^rvi 



+ (37 - 4)(1 - v'){l - S+) + (2 - 7)Sa^^"^ 



(5.30) 



+ [G+ - 2(7 - l)]Av, 
^apv"v^ = S+(i;2 + vl - 2vl) + VSE^ivl - vl) + 2V3E^V2V3 

+ 2y3(S3t;2 + S2t;3)t^i. (5.31 



Note that we have used the (Cc)i constraint ()2.52p to solve for r, giv- 
ing ()5.19|) . The (Cjs) constraint ()5.24|) is an immediate consequence of the 
definitions (12.1081) and = A/(3/32). Equation (jHIHI) is the /3-normahzed 
form of equation (j2.78|) . The evolution equations for Va are derived in the 
same manner as explained in Section 1^?^ for SH cosmologies. They were first 
derived by Weaver 2002 [SI] for generic G2 cosmologies. 
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Three useful auxiliary equations are 



g + 3S+ = 2(1 - A^) + diA - 8(^2' + S3') - |(2 - 7) 



(5.32) 



7-1 

7 



+ (37 - 4)f '(1 - S+) - i:apv''v'^ - 2(7 - l)v^Avi 
dtr = {q + 3S+)r + diq . 



(5.33) 
(5.34) 



Equation ()5.32|1 is derived from ()5.18p and ()5.21|1 , while ()5.34|1 is derived from 
the commutator [do, di] in ()2.118|) acting on /3. Equation ()5.32p becomes 
algebraic and convenient to use in the timelike area gauge ()5.3|) . Equation 
fl5.34|) is not used for numerical simulations, but is useful for analysis. 

It follows from (j2.66|) - (j2.7H) and (j2.117|) that the Weyl curvature variables 
are given by 



S+ = |S+ + |(E2' + E3') - |(S2 + S'J 

-lid,-r)A + liNl+Nl) + \^{v! 3 
^_ = i(l-3S+)S_--^(S3'-E2') 



+ lid, - r - 2A)N^ + |iV+iV_ + ::^^{vl - vl) 



— 1(1 — 3S_|_)Sx — ^23^2 



2^3 



1 -fn 



|(1 - V3S_)S3 - ^SxS2 - l:^viV2 



2G^ 

S2 = i(l + y3S_)S2 - ^SxSa - 



2G- 



-V1V3 



(5.35) 

(5.36) 

(5.37) 
(5.38) 
(5.39) 
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n+ = -Ar_S_ - A^xSx (5.40) 

= - r - A)Sx - |iV+S_ - S+Ar_ (5.41) 

= K^i - r - A)S_ - |iV+Sx - S+iVx (5.42) 

= -|(A + y3iVx)S2 - ^iV-S3 + ^^^3 (5.43) 

n, = liA - v^iVx)E3 + ^iV_S2 - ^^^2 . (5.44) 

Note that the {€0)2 and (^0)3 constraints have been used to simphfy H2 and 
H3 respectively. 

At this stage the /3-normahzed state vector is 

{E,\ S+, S_, Ex , S2, S3, N+, iV_, iVx , A R, n, V2, Vs, Qa). (5.45) 



The variable R is to be determined by choosing a spatial gauge, and the final 
form of the state vector depends on this choice, as we now describe. 

The spatial gauge 

We use a {t, x)-dependent rotation (see Appendix |B)) to rotate the spatial 
frame vectors {e^} so that S3 = 0. The evolution equation ()5.1(j|l and the 
constraint ()5.22|) then give 

E3 = 0, R= -v^SSx + (iV+ - V3N^)vi . (5.46) 

This is the shear spatial gauge ()4.64|) introduced for G'2-compatible SH cos- 
mologies. Alternatively, we can use the Killing spatial gauge ()4.63|) : 

= V3N_, R = -VsSx . (5.47) 

Both gauges specialize to ()4.65|) for generic G2 cosmologies with fewer than 
three tilt degrees of freedom. The motivation for using the shear spatial 
gauge rather than the Killing spatial gauge is the compatibility of the former 
with all G2-compatible SH cosmologies. 

Using the shear spatial gauge ()5.46|) . the final /5-normalized state vector 
for G2 cosmologies is 

X = (^/, S+, S_, Sx , S2, N+, N_,N^, A, n, vi,V2, Vs, Qa), (5.48) 
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together with the parameter 7. 

Separating Ei^ from the rest of the variables as follows, 



(5.49) 



where 



(5.50) 



we can write the system of evolution equations and constraints ()5.8|) - ()5.29p 
in a concise symbolic form. 
Evolution equations and constraints: 



In addition we write the expression (j5.19p for the spatial gradient r of the 
normalization factor f3 symbolically as 



Features of /5-normalization and separable area gauge 

The principal advantage of using /5-normalized variables and the associated 
separable area gauge is that it leads to evolution equations that are first or- 
der in the spatial derivatives as well as the time derivatives, i.e. a first order 
quasi-linear system of PDEs, which is desirable for numerical simulations. 
In contrast, if we were to use Hubble-normalized variables and the associ- 
ated separable volume gauge (see Uggla et al. 2003 |72j), then the evolution 
equations would contain second spatial derivative terms. The reason for this 
difference is that ()5.19p gives an algebraic expression for r, whereas when us- 
ing Hubble-normalized variables and the associated separable volume gauge, 
the expression for r contains first order spatial derivatives, which lead to 
second spatial derivatives since the equation for 9tS+ contains the spatial 
derivative dir. A minor variation of the formulation ()5.49p - ()5.54|) is possible 




(5.51) 
(5.52) 
(5.53) 




(5.54) 
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due to the fact that ()5.2H) is hnear in E+, allowing one to replace E4. by 



1-A^ + ^diA + A^n + 2A(A^xS_ - A^_Ex 
- E^ - E^^ - Es^ - - ^- - ^ - 



/m-A')] (5.55) 



(substitute (j5.19p in (|5.2ip and solve for E+). ^ Coupled with the fact that 
^iE_|_ does not appear in the evolution equations, this again leads to a first 
order quasi-linear system of PDEs. We shall make use of this possibility 
when doing numerical simulations. 



5.2 Invariant subsets of the /3- normalized state 
space 

In this section, we explain the important invariant sets and equilibrium points 
within the G2 state space, and the link with the SH state space. 

The G2 hierarchy 

We now revisit the hierarchy of the action of the G2 group and the tilt degrees 
of freedom, as shown in Figure and give the defining conditions in terms 
of the state variables ()5.50|) . 

For OT G2 cosmologies, ()3.5|1 and ()5.4fi|) imply 

E2 = 0. (5.56) 

For G2 cosmologies with one HO KVF, ()H.2j) and ()5.4f)|) imply 

= iV_ = = Ex . (5.57) 

For diagonal G2 cosmologies, ()5.5(i|l and ()5.57j) imply 

E2 = 0, iV+ = = = Ex . (5.58) 

The hierarchy of invariant sets is summarized in Table IKTT| which lists the 

^We restrict to models with < 1, i.e. models with a timelike gradient vector of the 
G2 area density. See van Elst et al. 2002 equation (62). 
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Table 5.1: The /3-normalized variables for G2 cosmologies. 

Invariant set State vector Y Additional restriction 

Generic G2 with 3 tilt degrees of freedom: none 

{T.j^,Y.^,T,^,T.2,N+,N^,Ny.,A,n,vi,V2,V3) 
Generic G2 with 2 tilt degrees of freedom: ^2 = iV+ = VSN^ 

S_, Sx , S2, A^-, iVx , A r;3) 
Generic G2 with 1 tilt degree of freedom: V2 = = = \/3A_ 

(S+,S_,Sx,S2,A_,A^x,A^^,^'i) 



OT G2 with 1 tilt degree of freedom: ^2 = ^3 = N+ = \/3A_, ^2 = 

(S+,E_,Ex,A^+,A_,Ax,A^^,^^i) 



One HO KVF with 2 tilt degrees of freedom: V2 = (A^+, A_, Sx) = 

(1;+, S2, Ax , A, 0,, vi,V3) 
One HO KVF with 1 tilt degree of freedom: V2 = V3 = (A+, A_, Sx) = 

(S+,S_,S2,Ax,A^^,i'i) 



Diagonal G2 with 1 tilt degree of freedom: V2 = V3 = (A+, A_, Sx, S2) = 
(S+,S_,Ax,A^^,^'i) 



variables and the additional restrictions. We have omitted in Table 15.11 
for simplicity and ease of comparison with Tables I4.1H4.4I for G'2-compatible 
SH cosmologies. 

SH dynamics again 

To specialize to G2-compatible SH models (see Figures IT^3.5|) . one sets 



d^Y = 0, r = 0. 



(5.59) 
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Then the separable area gauge coincides with the Gs-adapted gauge of 
SH cosmologies. The system ()5.52j) - ()5.54j) simplifies to 

dtY = g{Y) (5.60) 
C(Y,0) = 0, F(Y) = 0, (5.61) 

and (|5.51|) decouples. The system (|5.6(J|) - (j5.6ip . given in detail by (|5.8|) - 
()5.29|1 subject to ()5.59j) . describes the dynamics of G2-compatible SH cos- 
mologies, but in /5-normalized variables, and is fully equivalent to the Hubble- 
normalized system, given by ()4.26|) - ()4.46|) . 

The silent boundary 

The evolution equations ()5.8|) and ()5.34|) . 

dtEi' = iq + 3J:+)Ei\ dtr = {q + d,T.+)r + E^^d,q (5.62) 
imply that the conditions 

= 0, r = (5.63) 

define an invariant subset of the evolution equations, which we shall call the 
silent boundary. ^ 

It follows by setting Ei^ = and r = in ((^321)- (EMI), that the SH evo- 
lution equations ()5.60|) - ()5.61|) describe the dynamics in the silent boundary, 
but now with 

dxY unrestricted (subject to 1)5.6111 ). (5.64) 

In Section 15.31 we shall present evidence that the dynamics in the silent 
boundary approximates the dynamics of G2 cosmologies in the asymptotic 
regimes {t — > —00 or, A > and t — 00) in the sense that Ei^ and r tend 
to zero asymptotically. An important aspect of this approximation is that 
equilibrium points for the SH evolution equations are also equilibrium points 
for the G2 evolution equations on the silent boundary. 



^van Elst et al. 2002 "T^" define the silent boundary to he Ei^ = 0, with no restrictions 
on r. Using this definition, Andersson et al. 2004 [3] relate the vacuum G2 dynamics to 
the dynamics of spatially self-similar cosmologies (see sets (2) and (3) on page S39 of 
The difference here does not affect the past asymptotic dynamics, as spatially self-similar 
cosmologies and SH cosmologies have the same past asymptotic dynamics. 
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Figure 5.1: The Kasner parabola in the /3-normahzed state space. 
The Kasner equilibrium points 

On the silent boundary in the /?-normalized G2 state space, the Kasner so- 
lutions are described by a parabola of equilibrium points (see Figure EH)), 
given by 

n = o = nA, (iv+,iv_,iVx,ASx,S2) = , (5.65) 

and 

= 1 - 2E+ - . (5.66) 

The latter equation follows from ()5.2H) . From ()5.32j) . the deceleration pa- 
rameter q is given by 

q = -3E+ . (5.67) 
The parabola (j5.66|) is equivalent to the circle 

in the Hubble-normalized state space, as follows from ()5.73|1 below. The 
Taub Kasner point Ti (where /? = 0) is not represented on the parabola, 
because the /3-normalized variables are undefined there. The corresponding 
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transitions and Mixmaster dynamics can be described in the /^-normalized 
state space, but their representation in the Hubble-normahzed state space is 
simpler geometrically. 



The de Sitter equilibrium points 

On the silent boundary in the /^-normalized G2 state space, the de Sitter 
equilibrium points are given by 

nA = l, n = 0, = 0, N+ = N_ = = Ny, = A (5.69) 



or f ^ = 1 if 7 7^ I 

Vaix), 0<f2<l if7 = |. 



I (5.70) 



which are the same as (|4.i;i(i|) - (l4.i;i7|) . (|4.i;i9|) and (I4.142|) . except are 
functions of x rather than constants. From ()5.18|) . the deceleration parameter 
is given by 

q = -l. (5.71) 

The /3-normalized and the Hubble-normalized variables coincide here because 
S+ = (see equation (|5.72p below). 



Visualization of SH and G2 cosmologies 

In the SH state space, an SH cosmology at a fixed time is represented by a 
point, and its evolution in time by an orbit. A G2 cosmology at a fixed time 
can be represented by an arc of points in the SH state space, parametrized by 
the spatial variable x. The evolution of a G2 cosmology along one timeline 
is represented by an orbit in the SH state space, and the full evolution of a 
(j*2 cosmology is represented by a moving arc of points. We shall use this 
representation in Chapters IHl and IHl to illustrate G2 solutions in the SH state 
space. 

As discussed in Section 15.11 /3-normalized variables have the major ad- 
vantage of leading to evolution equations that are first order in the spatial 
derivatives, and will thus be used for analysis and numerical simulations. 
On the other hand, Hubble-normalized variables are directly related to the 
standard observational parameters, and give a simpler representation of the 
Kasner solutions. We shall thus use them for displaying the state space and 
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most of the plotted figures in Chapter IHl 

To avoid possible confusion between the notation of the Hubble-normalized 
variables in Section 14.21 and of the /3-normalized variables here, we shall use 
a superscript H for Hubble-normalized variables for the rest of Part |l] and 
Part HH of this thesis. Note that from 

Using this relation, /^-normalized variables can be converted into Hubble- 
normalized variables as follows: 

E?^^. «".^,etc. (5.73) 

In addition, the gradients of H and /3, defined by (^^-^i^ and 
are related according to 



rr r Ei^dxT,+ . 



It is important to note, however, that converting /3-normalized variables to 
Hubble-normalized variables does not change the temporal or spatial gauge 
(we are still using the separable area gauge and the shear spatial gauge). 

When the separable area gauge ()5.H) - ()5.2|) and the Gs-adapted gauge 
coincide, ^ the time variable t for the separable area gauge is related to the 
time variable for the Ga-adapted gauge according to 

= {\- . (5.76) 



5.3 Asymptotic silence 

There is considerable analytical and numerical evidence that the silent bound- 
ary plays a significant role in determining the asymptotic dynamics of G2 



■^This happens for SH cosmologies, and in the silent boundary of G2 cosmologies (see 
below). 
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cosmologies in the following sense: Ei^ and r tend to zero along almost all 
timelines as t ^ —00, and if there is a positive cosmological constant, as 
t — i> 00. The significance of this behaviour is that the asymptotic dynamics 
of G2 cosmologies will then be approximated locally, i.e. along individual 
timelines, by the dynamics in the silent boundary, i.e by SH dynamics. 

Past asymptotic silence 

We now define past asymptotic silence as 

lim = 0, lim r = 0. ^ (5.77) 

t—*—oo 

On the basis of numerical simulations and physical considerations about the 
existence of particle horizons (see below), we conjecture that generic G2 
cosmologies are asymptotically silent into the past: 

Conjecture 5.1 (Past asymptotic silence). A typical G2 cosmology is 
asymptotically silent into the past along almost all timelines. 

Some analytical evidence for asymptotic silence is provided by the evolu- 
tion equations. In the special case of vacuum OT G2 cosmologies with A = 0, 
the evolution equation ()5.8|) for Ei^ simplifies to 

dtE,' = 2E,' , (5.78) 

which gives 

= (Ei^)oe2* ^0 as t ^ -00, (5.79) 

where (i?i^)o is a positive function of x. We note that for a timeline whose 
orbit is past asymptotic to the Taub Kasner point Ti, t tends to a finite 
limit ts (see Section IHH]). So Ei^ tends to zero along orbits that are not past 
asymptotic to Ti. The evolution equation ()5.34|) for r also simplifies to 

dtr = 2r + E^^d^q . (5.80) 

^We note that, away from the Taub Kasner point Ti (where (3 — 0), Hubble-normahzed 
and /3-norniahzed variables only differ by a finite factor. Equation H5.77|l is thus equivalent 
to 

lim {Ei^)" = 0, lim r" = 0, 

since the limits are zero. 
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On the invariant set Ei^ = 0, this equation reduces to 

dtr = 2r, (5.81) 
and hence along orbits on the invariant set Ei^ = 0, r — ^ as t — —00. 

Future asymptotic silence 

Similarly, we define future asymptotic silence as 

lim Ei^ = 0, lim r = 0. (5.82) 

On the basis of numerical simulations, we expect that typical G2 cosmologies 
with A > are future asymptotic to the de Sitter solution. Evaluating (j5.8|) 
and (j5.34p on the de Sitter invariant set (j5.69p - (|5.7(J|) . we obtain 

dtE,' = -Er' , dtr = -r , (5.83) 

which imply that 

lim Ei^ = 0, lim r = 0. (5.84) 

Thus we expect that G2 cosmologies that are asymptotic to the de Sitter 
invariant set are asymptotically silent into the future. 

Physical interpretation of asymptotic silence 

The coordinate speed of light in the x-direction is given by 

^ = ±ArE,' , 5 (5.85) 



^The easiest way to see this is from the metric 

ds^ ^ -N^dt^ + {ei-)-'^dx^ + . 
For nuU geodesies in the x-direction, set ds — — dy — dz, and we obtain 

at 
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and it depends on the choice of temporal gauge and the parametrization 
of the t- and x-coordinates. In using the separable area gauge ()5.1|) - ()5.2|) . 
we have set J\f = 1. We also note that the slope of the light cone in the 
x-direction is the inverse of the coordinate speed of light: 

The physical speed of light has been set to unity by using geometrized units. 

If there exists a particle horizon in the x-direction, then its coordinate 
distance (from the observer at x = 0) at is given by 

XH{to) = I " Ei^Udt.^ (5.87) 

J —00 

The physical distance of the horizon at to from the observer is given by 

"'%(i„,.)£.'(<o,-)'"' 

where xnito) is given by (j5.87|) . 

With Af = 1, if a particle horizon exists, then x//(to) is finite and (as- 
suming a bounded derivative) Ei^ must tend to zero as t — > — cxd. If we also 
have r tending to zero, then we have past asymptotic silence. Thus, past 
asymptotic silence is interpreted as the existence of a particle horizon. 

Similarly, if there exists an event horizon in the x-direction, then its 
coordinate and physical distances are given by 

POO 

XEHito) = / Ei^Afdt (5.89) 

J to 
r-XEH(to) 

With A/" = 1, if an event horizon exists, then XEnito) is finite and (assuming 
a bounded derivative) Ei^ must tend to zero as t — 00. If we also have 



•^More precisely, XHito) is given by J_°^ Ei^{t, xughtit)) Af{t, xught{t)) dt, where xiight(i) 
is the path of the hght that reaches the observer at to. See Lim et al. 2004 ^Sj for the 
equations of the null geodesies. 
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r tending to zero, then we have future asymptotic silence. Thus, future 
asymptotic silence is interpreted as the existence of an event horizon. 

The significance of the spatial average of 1/Ei^ 

At a fixed time t, consider a spatial inhomogeneity with a coordinate wave- 
length of Ax = 2tt. Its physical wavelength is given by 

L(t)= / — — — -dx. (5.91) 

io /3it,x)E,^it,x) ^ ^ 

Suppose that a particle horizon exists. We can make a rough approximation 
of the ratio of this physical wavelength to the distance to the particle horizon. 
By assuming that /3Ei^ does not vary dramatically with x, it follows from 
(OHll and (linm that 

(5.92) 

dn xh 

Under a change of x-coordinate 

X = fix) , (5.93) 

Ei^ transforms as follows: 

E,' = f{x)E,' . (5.94) 
Note that the spatial average of (on x G [0, 27r] for example), given by 

is preserved under (I5.93j) provided that 

/(O) = 0, /(27r) = 27r, f'{x) > on s G [0, 2n]. (5.96) 

This average has physical significance, as we now explain. 
We choose /(x) such that 
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Assuming that Ei^ has an exponential growth rate from t = —00 up to time t 
(e.g. in the past asymptotic regime), ()5.87p reduces to xh ^ Ei^, and ()5.92p 
reduces to 




(5.98) 



The spatial average of is thus a measure of the ratio of the scale of spatial 
inhomogeneities and the distance to the particle horizon. The bigger Ei^ is, 
the shorter the wavelength is relative to the horizon distance. 

This understanding will be helpful in interpreting the numerical simula- 
tions in Chapter El in which we prescribe periodic boundary conditions by 
identifying x = 27r with x = 0. ^ In this situation the maximum coordinate 
wavelength is Ax = 27r. 



^In doing so, we do not regard the universe as being spatially finite, but rather as an 
infinite universe with repetitive patterns. In this way, we consider the particle horizon still 
exists even when xh > 27r. 



Chapter 6 

Explicit G2 cosmologies 



In this chapter we give two exphcit G2 solutions that can develop large spa- 
tial gradients on approach to the initial singularity. In the first solution, 
the asymptotic form of the (Hubble-normalized) state vector can have jump 
discontinuities or "steps", while in the second solution the asymptotic form 
can have "spike" discontinuities. The explicit solutions allow us to analyze 
the development of large spatial gradients in detail. 

6.1 An explicit solution with "steps" 

In this section, we discuss an explicit LRS G2 dust solution, which can de- 
velop step-like structures as it approaches the initial singularity. 

The framework 

We consider the class of LRS G2 cosmologies which are defined by the con- 
dition (jSSD: 

S2 = = S3 = S_ = Sx = iV„ = iV^, V2 = Q = (6.1) 
and we drop the index on t>i. We use the timelike area gauge (|5.3p : 

U = r, X=l, A = 0. (6.2) 
Recall that the spatial gauge is given by (|4.64|) . which now simplifies to 

N+ = 0, R = 0. (6.3) 

91 
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With these restrictions the evolution equations ()5.8j) - ()5.29j] simphfy dramat- 
ically. Firstly the Gauss constraint ()5.2H) defines S+ according to 



(6.4) 



and the remaining constraints are satisfied identically. Secondly, the expres- 
sions ()5.18|) and ()5.19p for q and r become 



1 + 3G;1[(7 - 1)(1 - v"") + ^v'']n - 3fiA 

1± (7- l)f2. 



"2G:"' 



The state vector ()5.48|) reduces to 

X={E,\Q,v,Qa). 



(6.5) 
(6.6) 

(6.7) 



Using (inill)-(inini) the evolution equations ^B>, (IK:^ and (IK:^ for 

these variables reduce to the following: 



9iE/ = (g + 3S+)Eii 



dtv + 



(7-] 



7 



i(l -3S+)(l+t;2 



1 - vydi \nn- [(37 - 4) - (7 - 1) (4 - 



7G'_ 

at(^A = 2(g+l)fiA . 



(2 - 7)G+r - 7t;[37 - 4 + 3(2 - 7)^+] 



(6.8) 



(6.9) 

-div 

(6.10) 
(6.11) 



An explicit solution 

Consider the ansatz v = 0, and with diQ 7^ 0. Equation ()6.6p implies r = 0, 
and then ()6.10|) implies 7 = 1 (otherwise Q = Q{t)). Equation ()6.5p simplifies 
to 

g=|(l-3fiA). (6.12) 
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Note that r = also implies (3 = P(t) and Qa = ^a(^), in accordance with 
equations ()2.108p and ()5.24p . Equations ()6.8p - ()6.15|) reduce to the following 
evolution equations: 

dtEi^ = {2-ln-3nA)Ei^ (6.13) 

dtn = 1(1 - n - 3nA)n (6.i4) 
dtnA = 3(1 - nA)nA . (6. is) 

These equations can be solved explicitly, yielding 

p-3t 

Q= -— , (6.17) 

(A£2 + e-3*) [1 + f{x)VAP + e-^'] 

' VA£2 + e-3* [1 + /(x)VA£2 + e"3*] 

where £ is a constant with dimension [length], and f{x) is an essential arbi- 
trary function, which is required to satisfy 



l + /(x)VA£2 > 0. (6.19) 

This condition ensures that > for t sufficiently large. If f{x) > 0, then t 
assumes all real values while if /(x) < 0, then t satisfies 

tsix) < t < oo, (6.20) 

where 

ts{x) = -iHfix)-' - Af) . (6.21) 

An arbitrary multiplicative function of x in ()6.18|) has been absorbed by 
means of x-reparametrization. The shear is given by 

^ VA£2 + [1 + /(x) VA£2 + e-3* ] ' ^ ' ^ 

If f[x) = const., the solution is an LRS Bianchi I dust solution with 
A > 0. In addition, if f{x) = 0, the solution is a fiat FL solution with 
fi + ^A = 1. 
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Structure of the singularity 

For simplicity we shall analyze the structure of the singularity using the 
explicit solution ()6.16|) - ()6.18|) . To show that an initial singularity exists, we 
show that the matter density fi 00. Using the definition of Q and and 
equations ()(i.l(jj) - ()(j.l7j) . it follows that 

^2^ ^ as t ^ if /(x) < 

^ 1 + f{x) y/XFT^^ [00 as t -00, if /(x) > 0. 

In order to determine the asymptotic state of the singularity, it is helpful 
to calculate the Hubble-normalized shear. Equations ()6.22|) and ()5.73|) lead 
to 



(6.24) 



It follows that 




for f{x) < (Taub Kasner) 
for f{x) = (fiat FL) (6.25) 
for /(x) > (LRS Kasner). 

Crucially, if f{x) has both signs, the limit of at the singularity has a 
jump discontinuity, i.e. a "step". Thus the solution can develop step-like 
structures as it approaches the initial singularity. 

One can gain insight into the creation of "steps" by thinking of the LRS 
G2 dust solution as an "inhomogenized" LRS Bianchi I dust solution. The 
Hubble-normalized state space for the LRS Bianchi I solutions is shown^ in 
Figure mm In terms of evolution into the past, the fiat FL equilibrium point 
is a saddle point, and the de Sitter equilibrium point is a source. The line 
going from the de Sitter equilibrium point to the fiat FL equilibrium point 
is the stable manifold of the fiat FL equilibrium point. At a given time, 
an LRS G2 dust solution is represented by a segment of points in the LRS 
Bianchi I state space. If f{x) > for all x, then > and hence the 
segment lies to the right of the stable manifold of the FL equilibrium point, 
and follows the arrows to the LRS Kasner equilibrium point as t ^ —00. 



^The state space follows from a routine analysis of equations lt).14|l - H6.15|l . and then 
making use of H6.4|l . 




Figure 6.1: The LRS Bianchi I state space (A > 0), and an LRS dust solu- 
tion as an "inhomogenized" solution at two different times. Arrows indicate 
evolution into the past. 
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Likewise, if f{x) < for all x, then the segment lies to the left of the stable 
manifold of the FL equilibrium point, and follows the arrows to the Taub 
Kasner equilibrium point as t ^ ts{x). Crucially, if f{x) has both signs, 
then the segment straddles the stable manifold of the FL equilibrium point. 
As t tends to ts{x) or —00, the segment is stretched out more and more. The 
part of the segment straddling the stable manifold of the flat FL equilibrium 
point is being stretched the most, giving a step-like structure in the spatial 
profile of S;^. 

This solution is a special case of the Szekeres solution (see Krasinski 
1997 im page 25], the case k = P' = 0). It is first given explicitly by Barrow 
& Stein-Schabes 1984 [3 case k = P' = 0]. They however did not discuss 
the past asymptotic dynamics. A special case of this solution, with A = 0, is 
first given explicitly by Bonnor & Tomimura 1976 12, Model PI with /i 7^ 0, 
P = 0]. They discussed the past asymptotic dynamics, but overlooked the 
possibihty that models can be past asymptotic to the Taub Kasner point. 

6.2 An explicit solution with "spikes" 

In this section, we obtain a new explicit vacuum OT G2 solution that develops 
"spikes" as t —y —00. Spikes are narrow spatial structures with steep slopes 
appearing on an otherwise mild profile of a variable of a cosmological solution 
(see Berger & Moncrief 1993 [H]). To generate this spiky solution, we apply a 
transformation introduced by Rendall & Weaver 2001 [02] to the Wainwright- 
Marshman solution (see Wainwright & Marshman 1979 [7^], Case I with 



The framework 

We consider the class of vacuum OT G2 cosmologies which are defined by 
the conditions 



(see ()5.56|) for the OT condition). We use the timelike area gauge ()5.3p : 



m = -^) 



n = 0, Va = 0, ^A = 0, S2 = 



(6.26) 



U = r, Af=l, A = 0. 



(6.27) 
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Recall that the spatial gauge is ()4.fj4j] : 

N+ = V3N_, R = - VsSx. (6.28) 

With these restrictions the evolution equations (j5.8|) - (j5.29j) simplify dramat- 
ically. Firstly the Gauss constraint 1)5.211) defines S_|_ according to 

and the remaining constraints are satisfied identically. Secondly, the expres- 
sions ()5.18|) and ()5.19p for q and r become 

g = 2 - 3E+ (6.30) 

r = -3(A^xS_ -iV_Sx) . (6.31) 

The state vector (|5.48p reduces to 

X=(Ei\S_,Sx,iV_,iVO. (6.32) 

Using (|6.26p - (|6.31|) the evolution equations (|5.8|) . (|5.12p - (|5.15|) for these vari- 
ables reduce to the following equations: 

dtEi^ = 2Ei^ (6.33) 

dtJ:- + Ei^d^N^ = 2^3(^2 - N^_) (6.34) 

+ = 2iVx (6.35) 

5tSx - Ei^d^N_ = -2v^(iV_iVx + S_Sx) (6.36) 

dtN^ - Ei^a^Sx = 2iV_ + 2y3(S_iV_ + A^x^x) . (6.37) 

Equation ()6.33j) can be solved explicitly to give 

Ei^ = e^^'^e'*, (6.38) 

and we reparametrize x to set g{x) = for simplicity. Since Ei^ ^ as 
t —00, the solutions are asymptotically silent (see Section 15. 2j) . The 
metric in this temporal and spatial gauge can be written in this form (see, 
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e.g. van Elst et al. 2002 jZH Appendix A.3]): 

ds^ = f3-'[-dt^ + {E,YV]+iy' [e''^''^\dy + Q{t,x)dzy + e-''^''^Uz^] , 

(6.39) 

where the variables iVx, Sx, N_) are defined in terms of P and Q (see 
van Elst et al. 2002 [7^ equations (180) and (181)], which are derived from 
the constraint (j2.8H) ): 

(S_,iVx) = ^^{dt,-E,'d,)P, (Ex, AT.) = -l^^e''{dt,E,'d,)Q. (6.40) 
The transformation method 

Rendall & Weaver 2001 [^2] introduce a two-step solution-generating method. 

The first step of the method is to perform the following transformation 
on the metric functions P and Q: ^ 

e^ = e-^[(Qe^)2 + l] (6.41) 

« = - W?FTT • ("2' 

The second step is to perform the so-called Gowdy-to-Ernst transforma- 
tion. Rendall & Weaver 2001 [^2]) define this transformation in terms of its 
effect on P and Q: 

P = -P -2t (6.43) 
Q, = e2(^+^*)Q. (6.44) 
4 = e'^Qi , (6.45) 

which generates a new solution,^ although in general Q cannot be obtained 
explicitly. 

We now describe the transformations in the framework of the orthonormal 
frame formalism. First, we show that the transformation ()6.4H) - ()6.42|) is 
a rotation of the frame vectors {e2, e^}, and describe it in terms of (3- 



^We have changed the sign of H6.42|l from equation (4) of |62| . 

^This is the same transformation as equation (7) of Rendall & Weaver 2001 ^ET, but 
looks slightly different because the time variable r in is related to our time variable t 

via r — —2t. 
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normalized variables. 

Consider a rotation of the frame vectors {e2, 63} (see Appendix |B|). 



(6.46) 



62 \ _ / coscp sm0 \ / e2 

63 / V -sin0 cos(j) J \ ea 

To preserve the spatial gauge condition ()6.28|) . it follows from ()B.2|) - ()B.7|) 
that the angle (f) must satisfy 

dt^= V3(-sin20S_ + cos20Sx) - VsSx (6.47) 
Ei^d^4> = V3(cos20 A^_ + sin20A^x) - VSN^ ■ (6.48) 

It turns out that these equations are satisfied if the rotation angle (p is related 
to P and Q according to 

as follows from differentiating ()6.49j) and using ()6.4()|1 . The transformation 
flB.6|) - ()B.7|) . with (p given by ()6.49|) . is in fact equivalent to performing the 
transformation ()6.4H) - ()6.42|) . This statement can be verified by differenti- 
ating ()6.4H1 - ()6.42|1 and expressing them in terms of (S„, A^x? Sx, A^-) using 
flOn|l . then comparing with (jR6|l - (jR7jl to obtain (jO^ . 

The Gowdy-to- Ernst transformation ()6.4Hj) - ()6.45j) is much simpler in terms 
of (S_, A^x,Sx, A^-), as follows from 

S_ = -S- - ^ (6.50) 

Nx = -iVx (6.51) 

Ex = A^- (6.52) 

iV_ = Ex . (6.53) 

It is easily verified that this transformation preserves the form of the evolution 
equations (Oijl - (I07I) . 

Composing the transformations ()B.6|) - ()B.7|) with ()6.49|) and ()6.50|) - ()6.53p 
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leads to 



S_ = —(cos 20 S_ + sin 20 Tix] 
iVx = sin 20 - cos 20 A^x 
Ex = cos 20 A^__ + sin 20 A^x 
= — sin 20 S_ + cos 20 Sx 



1 

V3 



(6.54) 
(6.55) 
(6.56) 
(6.57) 



where is given by ()6.49|) . 

In summary, given a solution (S_ 
ric functions P and Q are determined by ^6.40i )- Then (|6'.5^ )- |'6'.5?| ), with 
determined by P and Q according to \6.49j) , determines a new solution. 



,Arx,Sx,Ar_) of (fO^ - (f07D , the met- 



The Wainwright-Marshman solution 

Consider the ansatz 

Sx+A^_ = (iV-^0). (6.58) 

It follows from (jFnF)|l - (j07|l that S_ = iVx - Then imply 
iVx = 0. Solving (jOKjl then gives Sx = 756* /'(e^* - 2x), where / : R ^ M 
is freely specifiable, and e^* — 2x is the argument of the function /'. To 
recapitulate, we have 

^- = -^73' ^x=0' ^>c=-N- = j=/f'{e''-2x). (6.59) 

We can now obtain P and Q via ()6.40|) : 

P = -t, Q = f{e^' -2x). (6.60) 

This is the Wainwright-Marshman (WM) solution (see Wainwright & Marsh- 
man 1979 Hn], Case I with m = -^). ^ 



^The Gowdy-to-Ernst transformation leaves the WM solution essentially invariant - it 
changes only the sign of /. 
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The transformed Wainwright-Marshman solution 

We apply the transformation ()6.54|) - ()6.57|) to the WM solution and obtain a 
new explicit solution of the evolution equations ()6.33|) - ()6.37p : 

^- = -(-i73cos20+ ^sin20eV') - ^ (6.61) 

iVx = -^sin20e*/' (6.62) 

Sx = -^cos20e*/' (6.63) 

N- = ^ sin 20 + ^ cos 2(j) e'f (6.64) 

where 

cos 20=11^^^1^, sin20 = ^^^|^^, / = /(e^*-2x). (6.65) 
Structure of the singularity 

What solution(s) does the transformed WM solution tend to as t — >• —00? 
From (|6.65|) . we see that for a fixed x, 

lim cos 20 = < ^ /(-2a;) - 

^ \l iif{-2x)^0, ^ ' 

and sin 20 only pointwise because 

sin20 = ±l whenever /(e^* - 2x)e"* = ±1. (6.67) 

Then from ()6.61|) - ()6.65p . we obtain (assuming that /(e^* — 2x) is smooth): 

,~ ~ ~ ~ , f(-^, 0,0,0) for/(-2x) = 

' l(-i73'0'0'0) for/(-2x)^0. ^ ^ 



In terms of Hubble-normalized variables, we obtain 

;i,_4^) for /(-2a;) =0 
,^3, for /(-2a;) ^0 



:rrr. (6.69) 
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Figure 6.2: The two transition orbits shadowed by the transformed WM 
solution. (S:^,S:^) = (i,— ^y^) is marked by a square, and (S^,S:'^) = 
(y|, is marked by a diamond. 

i.e. the transformed WM solution tends pointwise to two distinct Kasner 
solutions if / has both signs. These two Kasner solutions are shown in the 
SH state space in Figure and are in fact connected by two A^:^ transition 
orbits. Recall that it follows from ()4.2H) that the two transition orbits 
lie on the sphere 

1 = (Ar^)2 + (s^)2 + (s;^)2^ (6.70) 

going from (S^, S^) = (f, to (tf, S^) = {%, as t ^ -00. We 

can think of the transformed WM solution as represented by a segment of 
points in the SH state space in Figure 16.21 If / has both signs, then this 
segment gets stretched along the two transition orbits as t ^ —00, thus 
creating spiky structures in the variables. 



CHAPTER 6. EXPLICIT G2 COSMOLOGIES 



103 




0<x<2ji 0<x<2ji 

Figure 6.3: Snapshots of the variables S_ and for /(•) = cos(«). 

To see this clearly, consider the example /(e^* — 2x) = cos(e^* — 2x). 
Figure inSl shows the snapshots of S_ and at t = 0, —1, —5. The variable 
S_ tends to for x = f , ^, ^, ^, and otherwise. The variable iV_ 

tends to zero pointwise, but not uniformly. The height of the spikes in N_ 
is We note that spikes also form in A^x and Ex, but these variables still 
tend to zero uniformly. 

In Section 18.31 we shall use both the WM solution and the transformed 
WM solution to test the accuracy of three numerical schemes. The trans- 
formed WM solution also provides a simple way to analyze the nature of 
spikes. We shall return to this in Section 
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6.3 Spiky and step-like structures 

Spikes and step-like structures are governed by the same mechanism. They 
have one thing in common - they form because the solution shadows two dis- 
joint branches of the unstable manifold (of an unstable equilibrium point). 
The meaning of a "trigger" , introduced in Section 14.11 can be extended to 
include the variable whose axis is parallel to the directions of the branches 
above. In the language of triggers, consider a trigger of an unstable equi- 
librium point. It assumes a certain value (we call it an "inactive value") at 
the unstable equilibrium point. Now consider a G2 solution in which the 
trigger of an unstable equilibrium point straddles its inactive value, e.g. the 
LRS dust solution in which straddles its inactive value of zero, and the 
transformed WM solution in which iV_ straddles its inactive value of zero. 
Then the trigger activates on both sides of the inactive value, but by conti- 
nuity there exists some point ^intercept in the middle where it cannot activate. 
Over time, a steep slope forms at xintercept- How the profile of the trigger 
extends further is determined by the unstable manifolds of the saddle point, 
whose orbits may lead directly to the sinks or indirectly through a sequence 
of saddle points. If these orbits are such that the trigger develops a monotone 
profile, then we see a step-like structure. Otherwise the profile resembles a 
non- monotone, spiky structure. Figure illustrates the initial smooth pro- 
file and the asymptotic formation of a steep slope at xintercept, and either a 
step-like structure or a generic spiky structure. 

We end this chapter with the following conjecture on what causes the 
formation of spiky or step-like structures. 

Conjecture 6.1. Consider a class of G2 cosmologies and their invariant 
subset of SH cosmologies. In an asymptotically silent regime, the G2 cos- 
mologies locally approach the attractor of their SH subset. Step-like or spiky 
structures will form if a G2 solution, represented by a segment of points in 
the SH state space, straddles the stable manifold of a saddle point. 
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Figure 6.4: Formation of spiky and step-like structures. 



Chapter 7 



Dynamics in the past 
asymptotic regime 



In this chapter, we make predictions about the dynamics in the past asymp- 
totic regime of G2 cosmologies, based on a heuristic analysis of the equations. 
Recall that we wrote the system (j5.8|) - (|5.29|) symbolically as (|5.51|) - (j5.54j) : 



where Y is the state vector in (j5.48|) . In analyzing the dynamics, we shall 
consider three types of behaviour: 

• the background dynamics, 

• the effects of saddle points, and 

• the effects of spatial derivative terms. 

The background dynamics is governed by the SH system ()5.60|) - ()5.61|) : 



which we regard as governing the dynamics in the silent boundary ()5.(jH|l : 



dtE^^ = B(Y)E,' 
dtY + M{Y)Ei'd,Y = g(Y) 
CiY,E,'d^Y) = , r = F(Y), 



(7.1) 
(7.2) 
(7.3) 



dtY = g(Y) 
C(Y,0) = 0, F(Y) = 0, 



(7.4) 
(7.5) 



El 







r = 0. 



(7.6) 
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We shall show that the saddle points of the SH system can combine with 
inhomogeneity to initiate some interesting effects in this system. Further- 
more, under the right conditions, the term M{Y)Ei^dxY becomes active 
and the system ()7.2|) exhibits some fascinating behaviour. We explain each 
type of dynamical behaviour in the sections below, and present numerical 
simulations that illustrate the behaviour in Chapter |H1 

7.1 Background dynamics: role of algebraic 
terms 

In Section we conjectured that a typical G2 cosmology is asymptotically 
silent into the past, suggesting that SH dynamics plays a major role in the 
past asymptotic dynamics of G2 cosmologies. The separable area gauge for 
G2 cosmologies also coincides with the G2 adapted gauge for SH cosmologies 
as t —>■ —00, since asymptotic silence entails r — > 0. 

Mixmaster dynamics 

With asymptotic silence, it is plausible that G2 cosmologies, like SH cosmolo- 
gies, also satisfy the BKL conjecture I and the Kasner Attractivity Conjec- 
ture, which will now refer to tilted Kasner circles )C±a as well as the standard 
Kasner circle /C. 

BKL conjecture I 

Along almost all local orbits of a typical G2 cosmology, 

lim n = 0, lim Qf, = 0. (7.7) 

Kasner Attractivity Conjecture 

Almost all local orbits of a typical G2 cosmology enter a neighbourhood of one 
of the Kasner circles. 

We emphasize that at this time there has been no progress in proving these 
conjectures for G2 cosmologies. They are, however, strongly supported by our 
numerical simulations. We believe that they are fundamentally important as 
regards the dynamics in the past asymptotic regime. 
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Once the G2 orbits enter a neighbourhood of one of the Kasner circles, 
Mixmaster dynamics takes over along each individual orbit, in accordance 
with SH dynamics. Proposition 14.11 concerning the occurrence of Mixmaster 
dynamics in G2-compatible SH cosmologies now leads directly to the follow- 
ing proposition. 

Proposition 7.1. Consider G2 cosmologies that satisfy the Kasner Attrac- 
tivity Conjecture. 

i) If the G2 is generic, Mixmaster dynamics occurs along typical time- 
lines} 

a) If the G2 is non-generic, Mixmaster dynamics does not occur. 

For generic G2 cosmologies, the background dynamics is that of SH cos- 
mologies of Bianchi type VI and VII. 

To recapitulate, the limits (j4.125|) - (j4.127j) for G'2-compatible SH cosmolo- 
gies should also hold locally (i.e. along typical timelines) for G2 cosmologies: 

Stable variables: lim {A,Ny.) = 0, (7-8) 

Non-overlapping triggers: \im {N22NSS, N22^x, NS3T.2) = 0, (7.9) 

t—>—oo 

Trigger variables: lim (A^_,E2,Sx) do not exist. (7.10) 

Again we emphasize that no progress has been made in proving Proposi- 
tion O and the limits (Tfini- dTTni) . The limits (I77|), (HH) and jlH), which 
assert that almost all local G2 orbits shadow the unstable manifolds of the 
Kasner equilibrium points, provide the foundation of Mixmaster dynamics 
in G2 cosmologies. 

Non-Mixmaster dynamics 

For more special G2 actions, there are too few triggers to sustain Mixmaster 
dynamics as some of the arcs on the Kasner circle become sinks (into the 
past) in the smaller state space. 

^We note that in contrast to the case of SH cosmologies, Nap = is not an invariant 
set in G2 cosmologies, due to the presence of the terms diY, and 5iSx in the dtNx and 
dtN- equations. 
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For OT G2 cosmologies, the background dynamics is that of OT SH 
cosmologies of type VI and VII. Thus we expect the local past attractor 
to be given by ()4.128|) : 

A- = }C aic{T^Qi). (7.11) 

By the phrase "local past attractor" we mean that almost all local orbits 
of a typical OT G2 cosmology are past asymptotic to the arc (jT.lip . As an 
example, the transformed WM solution tends to two distinct Kasner solu- 
tions, as given by ()6.69p . The local orbit of a spike point tends to the Kasner 
solution (S^, Y^^) = (y, — ^y^), which does not lie on /C arc(T3Qi), but other 

local orbits tend to the Kasner solution (E;^, S;^) = (y|, '~^^)-, which lies 
on /C aic{T^Qi). 

For G2 cosmologies with one HO KVF, the background dynamics is that 
of SH cosmologies of type VI with one HO KVF. We expect the local past 
attractor to be given by ()4.130|) : 

A- = [/C arc(S^ > -i(37 - 4)) n /C aIc{Q2T^QsT2) ] 

U [/C±i arc(S^ < -i(37-4)) n /C±i arc(g2Tig3T2) ] . (7.12) 

For diagonal G2 cosmologies, the background dynamics is that of diagonal 
SH cosmologies of type VI. Thus we expect the local past attractor to be given 

bv dorm : 

A' =IC arc(S:^ > -i(37 - 4)) U /C±i arc(S:^ < -i(37 - 4)). (7.13) 

Isenberg & Moncrief 1990 jH] proved that /C contains the past attractor 
for vacuum diagonal G2 cosmologies. Anguige 2000 iSj gave the asymptotic 
expansions for a generic class of perfect fluid diagonal G2 cosmologies that 
are past asymptotic to /C arc(S:^ > — 1(37 — 4)). 

7.2 Effects of spatial inhomogeneity I: role of 
saddles 

Recall that the state of a G2 solution at a fixed time is represented by an arc 
of points in the SH state space. We have seen in Chapter IHl that if this arc 
straddles the stable manifold of a saddle point then the solution will develop 
steep spatial gradients. In the language of triggers, if a trigger variable 
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Figure 7.1: The tilt variable v can develop step-like structures into the past. 

has both signs, then when the trigger activates, steep spatial gradients will 
develop in the neighbourhood of the point where the trigger variable is zero. 

In this section we discuss the steep spatial gradients in more detail, with 
a view to confirming the validity of BKL conjecture II. We shall sec that 
a mild initial inhomogeneity plus a zero in a trigger variable can lead to a 
strong asymptotic inhomogeneity, via local SH dynamics. 

Step-like structures 

The tilt variable v assumes a value close to or dzl when it is inactive, and 
undergoes transitions between and ±1 when it activates. If v has both 
signs on an arc of the Kasner circle K. when it activates, it develops step-like 
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structures, as illustrated schematically in Figure mi This is not unlike the 
step-like structures in the LRS dust solution in Section I6.H since the tilt 
transition will change t> = e>Otow = l — e, and f = — e < to f = — 1 + e, 
where < e -C 1. 

Spiky structures 

The curvature trigger and the frame triggers Sx and S2 assume a value 
close to zero when they are inactive, and become of order unity when they 
activate. If such a trigger has both signs on a Kasner arc when it activates, it 
develops spikes, as illustrated schematically in Figure ESI The spikes shown 
in Figure iniSl for the transformed WM solution are a simple prototype. 

In the case of the transformed WM solution ()6.61|) - ()6.65|) . the spikes are 
permanent - their amplitude does not tend to zero as t ^ —00. The spikes 
in the transformed WM solution form in a neighbourhood of isolated points 
where /(e^* — 2x) = 0, called "spike points", defined by 

/(e^* - 2a;spike(t)) = 0. (7.14) 

Since t appears in the argument of /, a spike point depends on t and moves 
toward a limiting spike point a;spike(— cxd) as t — > —00: 

a;spikc(t) = a;spikc(-oo) + |e^*. (7.15) 

Despite the fact that Ei^ —>■ exponentially fast, the steep slopes of the 
spikes ensure that M{Y)Ei^dxY in ()7.2j) is comparable to g{Y) in size, as 
shown by the asymptotic expansion of the transformed WM solution below. 

The decay rates as t — > —00 along the timehne x = a;spikc(— 00) are as 
follows (dropping tildes in (j6.6H) - (j6.64j) ): 

S_ = -^ + C(e2*) (7.16) 

= -2^/2 e^* + 0(e^*) (7.17) 

Sx = ^ /i e* + 0(e3*) (7.18) 

N_ = ^(/f - 1/2) e^* + 0(e5*) (7.19) 
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Figure 7.2: The triggers N_ and Sx can develop spikes into the past, as 
illustrated schematically in the above snapshots. 



CHAPTER 7. DYNAMICS IN THE PAST ASYMPTOTIC REGIME 113 



(7.20) 



(7.21) 

9tSx = ^/i e* + (9(6^*) ^lAT. = -^h e* + 0(e'*) (7.22) 

(7.23) 



where /i = /'(-2xspike(-oo)) and /a = /"(-2xspike(-oo))- 
For comparison, the following are the decay rates as t 
typical timeline x = a, where /(—2a) 7^ 0: 



V3_ 
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0(e- 
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_2^Ae2* + 0(e^*) 
A e* + (9(6^*) 
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A + ^ e* + 0(e^ 



3t\ 



-00 along a 

(7.24) 
(7.25) 
(7.26) 
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(7.28) 

e^* + 0(6^*) 
(7.29) 

g3t ^ (>^^g5t^) 

(7.30) 



= e^* + 0(e^*) , (7.31) 



where A = /'(-2a), A = /"(-2a) and /o = /(-2a). 

For typical points, from ()7.28j) - ()7.31|l . we see that the spatial derivative 
terms are negligible compared with the time derivative and algebraic terms in 
the evolution equations ()6.34|) - ()6.37|) . This supports the BKL conjecture II in 
Section IT!^ that the spatial derivatives in the EFEs for a typical cosmo logical 
model are not dynamically significant in the past asymptotic regime. 

For spike points, however, (|7.20p . (j7.22j) and (j7.23p show that despite 
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asymptotic silence, the spatial derivative terms are of the same order as the 
algebraic terms in the evolution equations ()6.34p . ()6.36p and ()6.37p respec- 
tively. This limits the validity of the BKL conjecture II within this class of 
models to typical timelines. 

For a general class of vacuum OT G2 solutions, Kichenassamy & Ken- 
dall 1998 equations (5)-(6)] gave an asymptotic expansion for the met- 
ric functions P and Q (see fl6.39p - ()6.40|) ) as t —00. Kendall & Weaver 
2001 |i62j transformed this expansion to produce expansions for P, Qt and 
Qx for a spiky solution. In order to obtain asymptotic expansions for the 
derivative terms in the evolution equations ()fi.34j) - ()(i.37j) . we found it neces- 
sary to increase the accuracy of their asymptotic expansions, as is done in 
Appendix O 

For typical points, the leading orders for the time and spatial derivatives 

are 



2 



dtt^ = e(e^'=* + e^^i-'^)*) diN^ = e(-te^*) (7.32) 

dtN^ = e{-te^') dit_ = e(e2*) (7.33) 

9tSx = e(e2(^-'=)*) diN. = 0(-te2(^+'')*) (7.34) 

dtN^ = e(e2'-*) 5iSx = 0(-te2('-'')*) , (7.35) 

where k{x) satisfies < A; < 1. We see that the spatial derivatives terms 
tend to zero as t ^ — 00 faster than the corresponding time derivatives in 
the evolution equations ()(i.34|l - ()fi.37|l . ^ 

For spike points, the leading orders for the time and spatial derivatives 

are 

dtt. = e(e^(^-'=)*) diNy, = e(e^(^^'=)*) (7.36) 

dtN^ = e(-te'*) dit^ = 6(6^*) (7.37) 

dtt^ = e{e^^^-'^') diN^ = e(e2(i-^')*) (7.38) 

dtN_ = e{e'"'' + e2(2-^0t) Q^f:^ = Q(_^g2(2-fc)t) ^ ^739) 

where k{x) satisfies < k < 1. Equations ()7.36p and ()7.38|) show that, despite 
asymptotic silence, the spatial derivative terms are of the same order as the 



■^The big 8 below is defined as follows: f{t) — Q{g{t)) means there exist constants ci 
and C2 such that f{t) satisfies cig{t) < f{t) < C2.g(t) for all t greater (or less) than to. 
•^^lE- has a marginal impact on dtNx- 
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algebraic terms in the evolution equations ()(i.34j] and respectively. 
On the other hand, ()7.39|) shows that diT^x dominates dtN^ in the case 
k > ^. This casts more doubts on the BKL conjecture II along the timelines 
corresponding to the spike points. 

It should be noted that the results (j7.32|) - (|7.39|) do not reduce exactly 
to the WM results because the WM solutions are special in a number of 
respects. For proper reduction, set k{x) = | in the more detailed expansions 
in Appendix O 

The main investigators of spikes (so-called Gowdy spikes in the context 
of vacuum G2 cosmologies) are Berger, Moncrief, Isenberg and their collab- 
orators. Berger & Moncrief 1993 ^T] first discovered spikes in numerical 
experiments. A more detailed numerical analysis is performed by Berger & 
Garfinkle 1998 HHI. 

7.3 Effects of spatial inhomogeneity II: role 
of spatial derivative terms 

In this section, we discuss inherently non-SH dynamics. 
Spike transitions 

Not all spikes are permanent. In Garfinkle & Weaver 2003 |^, numerical 
simulations show that for vacuum OT G2 cosmologies, spikes which develop 
during certain Kasner epochs are transient in the sense that they recede to 
a smooth profile via spike transitions. ^ Eventually, however, new spikes 
may appear at exactly the same spike point a;spike(i), and may be permanent. 
Unfortunately we have been unable to provide analytical heuristics to predict 
spike transitions. 

In Chapter ini we shall illustrate transient spikes in G2 cosmologies with 
generic G2 action, and describe spike transitions in more detail. We shall see 
that diN- is larger than at permanent spikes, enough to change the sign of 
9tSx and initiate a spike transition. 

^More precisely, spikes (as identified by their peak in the profile of S;?) that form near 

the Kasner arc spanning from T2 counter-clockwise to the point (S;^, S;^) = {—j, —^-7^) 
are transient. 
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Shock waves 

Shock waves can form if there are short-wavelength fluctuations with large 
enough amplitude in the perfect fluid variables. Heuristically, this can be 
explained as follows. For G2 cosmologies with one tilt degree of freedom, the 
evolution equation ()5.2(jj) for the tilt v has the following form: 



While other terms are not negligible, the left hand side of (j7.4()|l behaves like 
an inviscid Burgers' equation: 



It is well-known that solutions of Burgers' equation can develop shock waves, 
depending on the initial conditions (see, for example, Whitham 1974 |H3], in 
which Chapter 4 is dedicated to Burgers' equation). We predict that G2 
cosmologies can develop shock waves in the perfect fluid variables if the 
conditions are right, namely the occurrence of short- wavelength fluctuations 
with large enough amplitude in the perfect fluid variables {Q,v). In other 
words, the larger Ei^dxV is, the sooner v develops shock waves. 



dtv - 



[(37-4)-(7-l)(4-7H 



vEi dxV 



(7.40) 



dtu + udxU = 0. 



(7.41) 



Part II 

Numerical exploration of G2 
cosmologies 
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Chapter 8 

Equations for numerical 
simulations 

In the remaining two chapters, we show a variety of numerical simulations 
to support our analysis in Chapter [71 and to explore new phenomena. 

For simplicity we restrict our considerations to perfect fluid models with 
one tilt degree of freedom: 

V2 = V3 = 0, (8.1) 

although for some experiments we shall use vacuum models. In doing so, the 
shear spatial gauge ()5.46|) simplifies to 

S3 = 0, N+ = V3N_, R = -VsSx . (8.2) 

As discussed in Section 14.21 in the context of SH dynamics, the stabil- 
ity of the constraints is of great concern when doing numerical simulations. 
The concern is even greater for simulations of PDEs, which are less accu- 
rate than simulations of ODEs. Numerical experiments of G2 cosmologies 
show that the Gauss constraint 1)5.211) is the most unstable one during the 
transient regime, and the resulting error is unacceptably large. We remedy 
this problem by using ()5.2H1 to determine S+. The decision to use the (Cc)i 
constraint to determine r (see fl5.19|l ) also turns out to be equally important. 
The remaining (€0)3 and (C/j) constraints yield small, acceptable errors dur- 
ing the transient regime, and are stable in the asymptotic regimes. Ideally, 
the errors should be minimized, and the difficulties in doing so are discussed 
in Section lOl 

With the restrictions ()8.1|) - (|8.2p and solving (C^) = for S+, the state 



118 



CHAPTER 8. EQUATIONS FOR NUMERICAL SIMULATIONS 



119 



vector ()5.48j] reduces to the following for numerical purposes: 

X= (El^AS_,iVx,Sx,iV_,S2,fi,^;l,^^A), (8.3) 

together with the parameter 7, and we drop the index on vi. The full system 
of 10 evolution equations and 2 constraint equations is given in Appendix iDl 
We write the system of evolution equations ()D.3|) - ()D.11|) symbolically as 

5tX + M(X)9,X = (?(X). (8.4) 

Boundary conditions are required for simulations of PDEs. We shall pre- 
scribe periodic boundary condition, i.e. the numerical solutions are required 
to satisfy 

X(t, 0) = X(t, 27r) for all t. (8.5) 

It is customary to do so for numerical simulations of G2 cosmologies. 

By identifying x = 27i with x = 0, we have used up some but not all 
of the freedom in the choice of x-coordinate. A change of x-coordinate that 
satisfies ()5.9fj|l preserves this identification. We want to use all the freedom in 
the choice of x when specifying the initial conditions. Making Ei^ = const. 
in the initial conditions achieves this, and also makes the simulations more 
efficient. 



8.1 Specifying the initial conditions 

To perform a numerical simulation, we must first specify initial conditions 
that satisfy the boundary condition ()8.5p and the remaining constraints (Cc)3 
and (C/3): 

= (Cc)3 = ^/9,S2 - (r + - V3N^)J:2 (8.6) 
= (Cf,) = {E^d, - 2r)n^ . (8.7) 

A good way to specify the initial conditions is the following. We first 
specify the following eight variables: 



{E^\ r, A S_, AT,, Ex, N_, fi) 



(8.8) 
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on the interval x G [0,27r], and the parameter 7. ^ Because we are using 
the periodic boundary condition ()8.5|) . the initial data must be 27r-periodic. 
Then we calculate successively the following four variables 

(^A, E2, S+, Q), (8.9) 

where 

Q = ^n, G+ = l + {^-l)v\ (8.10) 
from the constraints ()8.7p . ()8.6|) . ()5.21|) and ()5.19p respectively: 

fiA = (^^A)oexp(^^ I^T^^) (S-ll) 

S2 = (S2)oexp[ I — dx j (8.12) 







E+ = I [1 - - - - f^fc - - I^a] (8.13) 
Q = |[-3AS+ -3A^xS_ +3iV„Sx -r] , (8.14) 

where {flA)o and (^2)0 are freely specifiable positive constants, and 

nk = Nl+Nl-l{Ei^d^-r)A + A\ (8.15) 

We often specify initial conditions which are so simple that Q\ and S2 can 
be computed by hand. We must make sure that ()8.1H) and ()8.12|) yield 
27r-periodic and E2, i.e. we must ensure that ()8.8p satisfies ^ 



2- 2r 



El 

2^ r + 3A - VSN: 



jdx = (8.16) 
dx = . (8.17) 



The tilt V is obtained by solving ()8.1()j) . which is a quadratic equation for v 



^Although during the simulation we use (|5.19|) to calculate r in terms of the state 
variables X, when specifying the initial conditions it is more convenient to specify r and 
use (|5.19|l to calculate Q and hence v. 

^If A = 0, then r does not have to satisfy H8.16|l . Likewise if S2 = 0, then r+3A— VSA^x 
does not have to satisfy H8.17|l . 
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giving 

2(Q/fi) 



We must make sure that ()8.8|1 and (T,2)o yield a v that satisfies 

-1 < < 1, or equivalently, - 1 < ^ < 1 . ^ (8.19) 

One way to ensure this is to specify a large Q, as we now explain. 

Consider the following format for the initial data, which we shall use in 
most numerical simulations in the next chapter: 

Ei^ = {Ei^)q , r = esinx, A = AQ + esmx, S_ = (S_)o + e sinx, 
TVx = v^Ao + sin X, = (iV_)o + esinx, Sx = (Sx)o + esinx, 

n = no + esmx, (8.20) 

where the subscript indicates a constant, and e is a small positive constant. 
Note that we have made the choice tq = and (A^x)o = V^^o in order 
to satisfy ()8.1(j|l - ()8.17|l . The remaining variables are determined by ()8.11|1 - 
(KT^i and (jOT^ . Equations (jKTT]) and (jO^ simplify to 

f^A = (^^A)oexp -2^^^cosx , S2 = (S2)oexp -Sj^fr^cosa; , (8.21) 



while S+ and v have complicated expressions. Expanding § as a power series 
in e, we obtain 

g = Ao + 2(iV_)o(Sx)o - 2v^Ao(S_)o 

- Ao [1 - AAl - {N^)l - (S_)^ - (Sx)^ - (S^)^] /n, + 0(e). (8.22) 

One can ensure f ^ < 1 by choosing a sufficiently small e and a sufficiently 
large ^q. It also helps to specify a small Aq (assuming Aq > 0). 



■^The second solution for v does not satisfy < 1. 

^One can prove this by plotting (|8.18|) on the box — 1 < g < 1, 1 < 7 < 2, and by 
plotting l|HTTn|l on the box -l<w<l, 1<7<2 using MAPLE. 
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Timelike area gauge A = 

If we use the timelike area gauge {A = 0), there is a simpler scheme. We first 
specify 

S_, iVx, Ex, iV_, n, v) (8.23) 
and 7. Then we calculate successively 

(r, ^A, S2, S+), (8.24) 

from (jOTI|l . (jHHI), dHini) and (jK:^^) respectively: 

r = -3(iVxS_ - iV_Sx) - l^t; (8.25) 
fiA = (^]A)oexp(^^ li^"^^) ^^-^^^ 



^ r 



E2 = (S2)oexp^ — ^^dxj (8.27) 
S+ = i [1 - - - - iV^ - _ ^ _ ^ ^g_28) 

where (f2A)o and (S2)o are freely specifiable positive constants. Note that 
in order to ensure that fl\ and S2 are periodic, equations ()8.16|) - ()8.17|) with 
A = must be satisfied. 



Vacuum subcase 

The scheme for the generic case does not apply to vacuum initial conditions, 
while the scheme for the A = subcase above does. We now provide a way 
to specify vacuum initial conditions with A > 0. We choose A^x = and 
specify 

iE,\ A, Sx, N_, S2). (8.29) 

Then the (Cc)3 and (Cq) constraints (j5.23p and (j5.2ip are solved for S+ and 
E_ in sequence, yielding 

S+ = l-(|Ei^9^1nS2-iV-Sx)M (8-30) 
S_ = [1 - - 2S+ - S2 - S2'] , (8.31) 
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where 

nk = Nl-l{Ei^d.^-r)A + A'^, r = -3ylS+ + 3A^_Sx . (8.32) 

8.2 Is the system hyperbolic? 

We ask the following question: does the system of equations (j8.4|) lead to a 
well-posed initial value problem? To answer this we rely on a theorem that 
states that if ()8.4|) is symmetric hyperbolic, then it leads to a well-posed 
IVP (see Friedrich & Rendall 2000 j24| pages 147-157]). By symmetric hy- 
perbolic we mean M(X) is symmetric, or more generally, that ()8.4|1 can be 
multiplied by a symmetric positive definite matrix A(X) so that y4(X)M(X) 
is symmetric (see Friedrich & Rendall 2000 [2H pages 148, 157]). One im- 
portant notion of hyperbolicity is strong hyperbolicity, which means M(X) 
is diagonalizable. Strongly hyperbolic systems are also symmetric hyperbolic 
in the general sense (see Friedrich & Rendall 2000 [211 page 161]). To our 
knowledge there are no other proofs of well-posedness for quasi-linear sys- 
tems of PDEs that satisfy other notions of hyperbolicity. See Gustafsson et 
al. 1995 i28i page 119] for these notions in the context of linear PDEs. 

The presence of dxA (contained in q and S+) in M(X) complicates the 
matrix. As a result we have been unable to determine whether the system 
()8.4j) is symmetric hyperbolic.^ We comment that if we use the timelike area 
gauge, where A is identically zero, then the system is symmetric hyperbolic,^ 
and thus has a well-posed IVP. Since we find that the numerical results 
obtained in this gauge and in the separable area gauge are qualitatively 
the same, we are confident that the separable area gauge is suitable for 
simulations. 



^The system is not strongly hyperbolic, because although M (X) has real eigenvalues, 
it is not diagonalizable on a set of points with measure zero. 

^This is proved by van Elst et al. 2002 [73] for OT G2 cosmologies. This result gener- 
alizes to the class of generic G2 cosmologies with one tilt degree of freedom, because the 
evolution equations for OT G2 are aumented by the dt'S'2 equation (jP.Sp . which has a zero 
characteristic speed of propagation. 
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8.3 Accuracy of numerical schemes 

In our numerical experiments, we make use of the following three numerical 
schemes: 

i) the fourth order Runge-Kutta method (RK4), 

ii) the 3-step iterative Crank- Nicholson method (ICN), and 

iii) the CLAWPACK code 

(see Gustafsson et al. 1995 |28| page 241], Teukolsky 2000 [70], and LeVeque 
2002 respectively). RK4 is a classic numerical solver for ODEs, but is 
suggested for application to PDEs by Gustafsson et al. ICN is a new method 
developed by Choptuik, but is not well-documented. CLAWPACK is developed 
by LeVeque, and is well-documented in 

We now explain the algorithm for each method briefly, and describe some 
tests of their accuracy when applied to the system ()8.4|1 . To conform with 
the notation in the literature on numerical algorithms for PDEs, we denote 
the state vector X by ?/ and rewrite the symbolic equation ()8.4|) as 

dty = f{y, d,y). (8.33) 

We denote the grid size by Ax and the step size by At. First, we have 
to decide how to discretize the spatial partial derivative d^y. The central 
differencing scheme for the derivative gives 

dMt, = — + 0{ (Ax)^), (8.34) 

which is second order accurate, i.e. the error is of order (Ax)^. The central 
differencing on five grid points is fourth order accurate: 



u X 4?/(t,x + Ax) -?/()!:, X- Ax) 

' ^ 3 2Ax 

1 y{t, X + 2Ax) - y{t, x - 2Ax) 



3 4Ax 



C( (Ax)^) (8.35) 



(see, e.g., Gustafsson et al. 1995 [2HI pages 19, 247]). 

Secondly, we have to give a scheme for time-stepping, i.e. to compute 
y{t + At, x), given y{t,x). Let ?/„ and yn+i denote the computed values for 
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y{t) and y(t + At) respectively. The fourth order Runge-Kutta method (RK4) 
is described as follows (see, e.g., Gustafsson et al. 1995 |2Hl page 241]). For 
brevity we write dxUn) as /(?/„) (or think of / as an operator). 

ki = AtfiVn) 

k2 = AtfiVn + ifci) 

ks = AtfiVn + \k2) 

ki = Atfivn + ks) 

Vn+l =yn+\{kl + 2k2 + 2^3 + ^4). 

Using Taylor polynomials, it is straightforward to show that RK4 is fourth 
order accurate (the error in the spatial derivatives may for the moment be 
ignored) : 

Vn+i = y{t + At, x) + 0{ {Atf), (8.36) 

i.e. the power of At in the error term is four orders higher than that in 
the formula for Un+i- Now combining RK4 with the 5-grid central difference 
scheme, one obtains a fourth order accurate finite difference scheme. 

The 3-step iterative Crank-Nicholson method (ICN) is described as fol- 
lows: 

ki = Atfiyn) 

k2 = Atf{y^ + ifci) 

ks = Atfivn + 1^2) 

Vn+l = Vn + ks. 

ICN is second order accurate.^ Now combining ICN with the 3-grid central dif- 
ference scheme, one obtains a second order accurate finite difference scheme. 
Since RK4 evaluates the equations four times at each time step while ICN only 
three, it takes RK4 4/3 times as long to run as it takes ICN. 



^See Teukolsky 2000 While an n-th order Runge-Kutta method is n-th order 

accurate, Teukolsky's analysis shows that all ICN methods are only second order accurate, 
and the 3-step ICN method is the cheapest stable method among them (see Teukolsky 
2000 ,70: )- We note that there are more sophisticated implementations of the iterative 
Crank-Nicholson method (Choptuik 2004 2S|). Teukolsky claims that there is no benefit 
in applying more iterations. Choptuik, however, has reservations about this claim, on the 
basis that Teukolsky's idealized analysis is incomplete. 
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The CLAWPACK code (see LeVeque 2002 [ISj) handles the evolution equa- 
tions ()8.4|) in two stages for each time step, using the so-called Godunov 
splitting (see LeVeque 2002 Pl pages 380-388]). In the first stage, CLAWPACK 
solves the following system (referring to ()8.4|) ): 

dty + M{y)d,y = 0, (8.37) 

where M{y) must be diagonalizable. i.e. it can be written as 

M = QDQ~\ (8.38) 

where Z) is a diagonal matrix of real eigenvalues, and Q is a matrix con- 
structed from the corresponding eigenvectors. The upwind finite volume 
method (see LeVeque 2002 page 73]) is then used to discretize (|8.37p . 
To explain the discretization, consider rewriting ()8.37p as 

dtw + Dd^w = 0, w = Q-^y, (8.39) 

pretending that Q is constant. Since D is diagonal, this system decouples 
into independent advection equations, and the eigenvalues determine the 
velocities (and hence the upwind directions) of the waves in the characteristic 
eigenfields w. So we are able to determine the upwind direction for each wave. 
In finite volume methods, an x-interval is divided into cells, and fluxes are 
computed at the interface of two adjacent cells. To obtain the data at the 
interface, some averaging scheme must be used (see LeVeque 2002 03 Section 
15.3]), and here we simply take the arithmetic average. The discretization of 
(Pn?7|) yields 



iv: -y?-i} + iy?+i - y?) , (8-4o) 



where 

Af\=Ql,Df\{Ql,)-\ (8.41) 

*~2 '2*2 2 

and D = + , where contains only positive eigenvalues while D~ 

contains only negative eigenvalues. The vector A'^_i{yi — yi-i) / Ax describes 

* 2 

fluxes to the right from the {i — l)-th cell to the i-th cell, etc. For details see 
LeVeque 2002 [13 pages 78-82]. The first stage is first order accurate. 
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In the second stage, CLAWPACK solves the system 

dty = g{y). (8.42) 

This part can be discretized by any ODE solver. We solve this part by using 
the explicit Euler method, 

Vn+l =yn + ^tg{yn), 

which is first order accurate. 

Overall, CLAWPACK is first order accurate. While RK4 and ICN do poorly 
in simulating shock waves (Taylor polynomials do not apply at points of 
discontinuity), CLAWPACK sacrifices high accuracy for smooth solutions for 
the ability to simulate shock waves reasonably well. The error analysis of 
CLAWPACK near shock waves is based on conservative systems of the form: 

dty + d,iGiy)) = giy). 

However, the system ()D.3|) - ()D.11|) is not conservative, due to the presence of 
inhomogeneity in Ei^ and the non-conservative pair (In ^2, arctanh v). This 
will likely reduce the accuracy of CLAWPACK when simulating shock waves. 
Another more serious abuse of CLAWPACK is that we have put dxA (contained 
in q and E+) in the right hand side of the equations {g{y) becomes g{y, dxy)),^ 
which makes it unsuitable for shock waves simulations. The best we can do 
is to simulate shock waves using the timelike area gauge A = 0. The abuse 
becomes relatively harmless if Ei^ is small. For example, CLAWPACK does a 
better job than the other two schemes at simulating permanent spikes under 
low resolution (see Section IHTij) . so there is value in keeping CLAWPACK. 

To measure the rate of convergence for each method, we make three runs, 
using 512, 1024 and 4096 spatial grid points. The data using 4096 spatial 
grid points serves as the reference, and the errors of the two runs with cruder 



^CLAWPACK requires a coefficient matrix M(y) that is diagonalizable. By putting dxA 
in the left hand side, CLAWPACK is liable to crash due to unbounded determinant of the Q 
matrix in H8.88|l . 
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grids are compared. We use both L2 and Loo norms to measure the errors. 



Table 8.1: The rate of convergence of numerical methods: errors relative to 
the solution on 4096 grid points as measured in the L^o norm. 



Method 


error (512) 


error (1024) 


ratio 


convergence 


RK4 


2.9746 X 10"^ 


1.9852 X 10-^ 


14.9839 


fourth order 


ICN 


2.2908 X 10-3 


5.3947 X 10-^ 


4.2464 


second order 


CLAWPACK 


1.7974 X 10-^ 


7.8783 X 10-2 


2.2814 


first order 



In Table 18.11 we present the errors in L^o norm. A method is said to be 
n-th order accurate in time if its error (approximately) shrinks by a factor 
of 2-" when the number of grid points is doubled. Errors in L2 norm are 
of similar size and are not presented. The initial data is chosen to be (|8.20|) 
with 7 = 1 and 

e = 0.01, {Ei')o = l, Ao = 0.2, (S_)o = 0.2, (Sx)o = 0.4, 
(iV_)o = -0.4, fio = 0.5, (f]A)o = 0.7, (S2)o = 0.1. (8.43) 

The simulations run from t = to t = —1. Longer simulations will see 
the formation of spikes. The order of accuracy is derived from Taylor poly- 
nomial approximation. With decreasing resolution around the spikes, the 
approximation becomes poorer and poorer, and eventually becomes totally 
useless. Higher resolution will be needed to prolong simulations. The CFL 
stability condition (see Gustafsson et al. 1995 |28» page 54]) is set at At < 

^Let be the m-th component of the vector y evaluated at the i-ih grid point. Then 
the L2 norm of the error in y is given by 

lly - yrcf II2 = ./^E(yr-(yrcf)r)' , 

y z,m 

where N is the number of grid points. The Loo norm of the error in y is given by 

\\y - yrcf lloo = max I?/™ - (2;rof)ri • 
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Ax/ maxEi^. This serves as a cap for the step size At to maintain the stabil- 
ity of the solution. The preferred value for At is set at At = 0.9Ax/ maxE'i^, 
i.e. the step size is computed as follows: 

At = 0.9Ax/maxEi^(t,x). (8.44) 

After the data at time t + At is obtained, we check that 

At < Ax/maxEi^(t + At,x). (8.45) 

If this is satisfied, we continue; otherwise we recalculate the data with a 
smaller step size 

Atsmaiicr = 0.9 Ax/ max ^ (t + At, x) . (8.46) 
This smaller step size will satisfy the CFL condition: 

Atsmallcr < Ax/ maxEi^(t + Atsmaller, x) . (8.47) 

Now since the reference solution is computed by the same method, we can 
only conclude that each method converges at the respective rate. Do they 
all converge to the real solution? To answer this we need an exact solution 
against which the error of the methods can be measured. Unfortunately 
there are no exact generic G2 solutions available thus far. The most suitable 
possible test solution is the WM solution (see Section lU?^ . which is an exact 
vacuum OT G2 solution, containing one arbitrary function / : M ^ R. It is 
given by (!(i.59jl: 

S_ = -^, N^=0 = A, Sx = -iV_ = ^eV(e'*-2x), E,' = e^K 

(8.48) 

For the test, we choose 

/(e^* - 2x) = sin(e2* - 2x) with < x < 27r. (8.49) 

The solution then satisfies periodic boundary conditions. The simulations 
run from t = Otot = — 1. The CFL condition is set at At < Ax/ maxEi^, 
and the preferred value for At at At = 0.9Ax/ maxii^i^. In Table IH^ we 



'We also impose an absolute cap to the time step: At < 0.1 for stability. 
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present the errors in Ex in the L^o norm. Errors in the L2 norm are similar 
and are not presented. 

Table 8.2: The accuracy of numerical methods: errors in the simulated Sx 
for the WM solution as measured in the L^o norm. 

Method error (512) error (1024) ratio accuracy 



RK4 
ICN 
CLAWPACK 



9.5834 X 10-9 
9.8275 X 10-*^ 
4.3212 X 10-3 



6.2809 X 10-^^ 
2.7601 X 10-6 
2.2423 X 10-3 



15.2581 
3.5605 
1.9272 



fourth order 
second order 
first order 



Table 8.3: The accuracy of numerical methods: errors in the simulated E_ 
for the transformed WM solution as measured in the L^o norm. 



Method 


error (512) 


error (1024) 


ratio 


accuracy 


RK4 


1.5157 X 10-5 


9.7292 X 10-^ 


15.5789 


fourth order 


ICN 


2.6242 X 10-3 


6.6930 X 10-^ 


3.9208 


second order 


CLAWPACK 


2.4204 X 10-2 


1.2511 X 10-2 


1.9347 


first order 



The WM solution can be transformed into a new solution (see Sec- 
tion 123). We test the codes against this solution, using the initial condition 
for the WM solution above. In Table IH751 we present the errors in S_ in L^o 
norm. 

The tests confirm that the methods produce solutions that converge to 
the exact solution at the expected rates, which gives us some confidence that 
the methods are reliable. 
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8.4 Unsatisfactory aspects of the numerical 
simulations 

Instability of (Cc)3 

The evolution equation for {€0)3 is given by 

dt{Cc)3 = (2g - 2 + 3S+ + y3S_)(Cc)3. (8.50) 

This constraint is stable into the past near Kasner equilibrium points. How- 
ever, when S2 is large, (Cc)3 grows as t ^ —00. This is a cause for concern 
especially when {€0)3 grows at the beginning of the simulation (i.e. before 
entering the asymptotic regime). Using finer grids or taking smaller time 
steps reduces the size of (Cc)3, but is very inefficient. Unfortunately, we find 
that adding the constraint in the evolution equations does not change its 
instability. Constraint instability is a recurring problem in numerical rela- 
tivity. For other methods in controlling constraint instability, see Lindblom 
et al. 2004 jHO]. 

Resolution 

As we have not implemented an adaptive mesh refinement scheme, high- 
resolution simulations are expensive for the numerical schemes we have. This 
imposes a severe limitation when spikes occur. Without mesh refinement, we 
can only adequately simulate spikes for a short time before the spikes become 
too narrow to resolve numerically. If Ei^ is sufficiently large near a spike, RK4 
can generate spurious oscillations, while ICN also generate some overshoot. 
Only CLAWPACK handles spikes reasonably well in all situations. 

The dust case 

For dust (7 = 1) the equations are not strongly hyperbolic (the matrix M(X) 
is not diagonalizable at any point in spacetime), regardless of whether A = 
or not. We shall not simulate the dust case in this thesis. 



For example, 

dt^2 = (9 - 2 + %/3S_)I]2 + (Cc)3 • 
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The Taub Kasner point Ti 

The Taub Kasner point Ti given hy P — causes difficulties for numerical 
simulations. It causes blow-up in the /?- normalized variables because f3 is zero 
here. This problem cannot be avoided by using Hubble-normalized variables 
or other gauges, because the Taub Kasner point is qualitatively different from 
other Kasner points, namely it is non-silent in the x-direction. We can only 
simulate solutions that do not get too close to the Taub Kasner point during 
the period of interest. As a general rule of thumb, the more inhomogcneous 
the initial condition is, the more likely we will run into the Taub Kasner 
point at an early stage in numerical simulations. 



Chapter 9 

Numerical simulations 



In this chapter, we present a variety of numerical simulations of G2 cos- 
mologies that confirm the analytical predictions made in Chapter [7| The 
numerical results also illustrate the spike transitions, which we were unable 
to predict analytically in Chapter [7| In addition we describe simulations of 
the future asymptotic approach to de Sitter, and the close-to-FL regime. 

For most of the simulations, ^ we shall specify initial conditions by fol- 
lowing the format (j8.2m) introduced in Chapter |S| 



El" = {Ei^ 



) 



VsAq + ^ sin X, 



esinx, A = Ao + esmx, E_ 
= (A^_)o + esina;, Sx 
1] = i7o + e sin X, 



[S_)o + esinx, 
(Sx)o + esinx, 

(9.1) 



where the subscript indicates a constant. and S2 are given by ()8.21|) : 



-2 



(£^1^)0 



cos X 



S2 = (S2)oexp 



-3 



cosx 



(9.2) 



while S+ and v are computed numerically using ()D.13|) . ()8.14p and ()8.18|) . 
and have complicated expressions. For more special G2 actions, we set the 
appropriate variables to zero. 

While /5-normalized variables are used in the simulations, Hubble-normalized 
variables are displayed in most of the plots. 



^The exceptions are the simulation of vacuum models in Section [9.2.3l the shock wave 
simulation in Section l9.4.2l and the close-to-FL simulation in Section 
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9.1 Asymptotic silence 

In this section we shall provide evidence to support the occurrence of asymp- 
totic silence not only as t — — oo, namely, 

lim Ei^ = , lim r = , (9.3) 

t— »— oo 4— »— oo 

but also as t ^ oo in the presence of a cosmological constant. 

Using the format (jfj.lll for initial conditions, consider the following initial 
data at t = 0, with 7 = |: 

e = 0.01, (^i^)o = l, Ao = 0.2, (S_)o = 0.2, (Sx)o = 0.4, 

(iV_)o = 0.4, fio = 0.5, (fiA)o = 0.7, (S2)o = 0.1. (9.4) 

We ran ICN with 512 grid points, from t = to t = —50, and from 
t = to t = 5 storing data at intervals of t = 0.1. First, we plot the 
maximum and minimum of In(ii^i^)^ and In |r| against t in Figure I^TTl We 
see that (£"1^)^ and r decay exponentially both into the past and into the 
future, supporting the conjecture that the past and future dynamics are 
asymptotically silent.^ Thus it is reasonable to say that the local background 
dynamics of G2 cosmologies is the SH dynamics. 

9.2 Background dynamics 

In this section we illustrate the role of SH dynamics in G2 dynamics. 
9.2.1 Local Mixmaster dynamics 

Still using the initial data ()9.4j) . we plot the maximum and minimum of 
ln(f2^) and ln{Q^) against t in Figure IHUl We see that and decay 
exponentially into the past, supporting the BKL conjecture I. 

Fi gures f9 . 3H9 . 41 show the evolution along a fixed timeline x = const, (the 
grid point 100 of 512). They show the generalized Mixmaster dynamics into 
the past and the approach to de Sitter into the future. 



^Actually, this is only true along almost all timelines. In Section [9.4.11 we shall see 
that along some timelines r can become large, but the spatial resolution in this run is too 
low to simulate it correctly. 
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Figure 9.1: Exponential decay of {Ei^)^ and r as t —oo and as t ^ oo, 
showing the maximum and minimum values of \n{Ei^)^ and In |r | . The graph 
for min In Irl is broken because r has both signs. 

Q<x<2'K 
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t=o 

* t=5 




Figure 9.3: Typical evolution along a fixed timeline, showing the transition 
orbits. 

Evolving into the past from t = 0, the solution enters the Kasner regime 
at about t = —3. One prominent S2 transition is observed at about t = 
—9, followed by an alternating sequence of and Sx transitions (see also 
Figure 19.41 below) . 

Figure 19.41 shows that v is stuck at v = last^— 00, and tends to a 
no n- zero limit as t — > 00. For more details about the future asymptotics, see 
Section 1^31 

Table inm shows the expected and actual growth rates into the past (posi- 
tive value indicates growth into the past) along grid point 100, at two selected 
times. The actual rates are in good agreement with the expected rates, which 
can be inferred from the /3-normalized version of ()4.89|) - ()4.91|) . 
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Figure 9.4: Typical evolution along a fixed timeline. 
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Table 9.1: Predicted and actual growth rates of the triggers into the past 
(positive value indicates growth into the past) along grid point 100 of 512. 



Triggers 


Sx S2 VI -v^ 


At t = -50 (0 = 1.2076 or 69.2°) 


Expected rate 
Actual rate 


-4.3894 2.3894 -0.4083 0.7864 
-4.4174 2.3920 -0.4083 0.7864 


At t = -35 (0 = 4.1341 or 236.9°) 


Expected rate 
Actual rate 


4.3976 -6.3976 -0.4173 -3.6161 
4.4040 -6.4323 -0.4173 -3.6161 



9.2.2 Spatial variation in the background dynamics 

The timing of the transitions in Mixmaster dynamics varies from one timeline 
to another. Continuing with the initial data ()9.4|1 . Figure PTHl shows some 
snapshots of the triggers and , illustrating this variation in the timing. 
The variation in the timing leads to a pulse-like profile of the triggers. 

As discussed in Section 15.31 the particle horizon is given by Ei^dt. 

Since ii^i^ = 1 at t = here and decays like e^*, Ei^dt is approximately 
equal to |, and so information can travel only a finite distance of Ax ~ | 
as t —>■ 00. But the pulses in Figure 1^31 travel more than Ax = |, which is 
only possible if they propagate faster than light. How is this possible? In 
fact, the apparent superluminal speed is an illusion created by the variation 
in the background SH dynamics from one timeline to another (e.g. variation 
in the timing of the transitions in Mixmaster dynamics). An analogy to this 
phenomenon is a Mexican wave in a stadium. 
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t=-10.8 




0<x<2jt 0<x<27t 

Figure 9.5: Snapshots of triggers and , showing pulses apparently 
travelling at superluminal speeds. 
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9.2.3 Multiple transitions 

It is simplest to analyze multiple transitions in the absence of spikes or step- 
like structures. Since it is difficult to produce a spike-free solution over a 
very long run in a perfect fluid solution,^ we resort to simulating vacuum 
solutions. 

Using the scheme for the vacuum case in Section 18.11 we choose = 
and specify the following vacuum initial condition: 

Ei^ =A = 0.1 + 0.01 sin 
Ex = -A^- = 0.8 + O.lsina;, S2 = 0.1 . (9.5) 

Then E+, S_ and r are given by ((1301) -(IH32)- 

We ran CLAWPACK with 512 grid points, from t = to t = —50, storing 
data at intervals of t = 0.1. Figure shows bands of level curves indicating 
the activation of the triggers S2, Sx and N_. The result is spike-free up to 
t = —33. Spikes occur in the lower right quarter of Figure I^H)! 

In Figure 19.51 of the previous section, we saw that there is variation in 
the timing of the transitions in Mixmaster dynamics, due to spatial inhomo- 
geneity. Figure 19.61 shows that the sequence of the transitions can also vary 
from one timeline to another. For example, proceeding along the timeline 
X = 71 as t —00, the solution undergoes an transition ffist at about 
t = —10 before undergoing a S2 transition at about t = —20. In comparison, 
proceeding along the timeline x = ^ as t — — 00, the solution undergoes a 
S2 transition ffist at about t = —17 before undergoing an iV_ transition at 
about t = —18. Thus, by continuity, there must be timelines between x = tt 
and X = ^ that undergo both E2 and N_ transitions at the same time, i.e. 
they undergo a multiple transition. 

The intersection of the bands indicates multiple transitions. We see that 
such intersections are inevitable in the global picture. Thus multiple tran- 
sitions will occur from time to time globally, as t —>■ —00. Locally along a 
fixed timeline, it is possible that the probability of multiple transitions occur- 
ring will tend to zero as t — > —00, in agreement with the multiple transitions 
conjecture in SH dynamics, but we cannot make a strong claim in this regard. 



■^That the initial condition for fl is large (see the comment following H8.19|l 'l makes it 
difficult to predict whether a solution will be spike-free and will stay away from the Taub 
Kasner point for a long period of time. 
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< X < 2jt 



Figure 9.6: The level curves of = 0.1 (solid), = —0.1 (dashed) 
and = 0.1 (dotted), indicating the activation of the triggers. Multiple 
transitions occur when two bands overlap. For example, in a neighbourhood 
of t = —32 and x = 7r/4, the solid band and the dotted band overlap, 
indicating a multiple transition involving and S^. 
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9.3 Effects of spatial inhomogeneity I: role of 
saddles 

In this section, we illustrate the role of saddle points in creating step-like 
or spiky structures. We also take this opportunity to illustrate the past 
attractors of diagonal and OT G2 cosmologies. 

9.3.1 Step-like structures 

We now illustrate step-like structures in the tilt variable v in diagonal G2 cos- 
mologies, and provide evidence to support the claim that the past attractor 
is given by ()7.13|) : 

A- = }C arc(S:^ > -i(37 - 4)) U /C±i arc(S:^ < -i(37 - 4)). (9.6) 

Recall that diagonal G2 cosmologies satisfy ()5.58|) : 

S2 = = Sx = TV- . (9.7) 

Consider the following initial data at t = with 7 = | (in conjunction 
with (jnH) and (ETj)): 

e = 0.3, (^/)o = l, Ao = 0, (S_)o = 0.2, Ex = 0, 

N_ = 0, fio = 0.6, {nA)o = 0.1, S2 = 0. (9.8) 

We ran CLAWPACK with 512 grid points, from t = to t = — 10, storing 
data at intervals oft = 0.1. Figure I^TI shows the snapshots of the tilt variable 
V at t = —4, t = —5 and t = —10. We see the strong suggestion that three 
step-like structures develop as t — » —00: one between v = ±1, and two 
between v = and v = 1. The one between v = ±1 involves an unstable 
state V = 0, while the two between v = and v = 1 involve a string of Kasner 
equilibrium points with E;^ = and — 1 < f < 1 that connect the stable arcs 
JC arc(S^ > 0) and /C±i arc(S:^ < 0). 

Figures 1^^9.101 show the snapshots of the variables (E^,S:^,f) in the 
(E;^, T,^, v) space. Each dot represents the state at a grid point. We see that 
the state vector approaches the past attractor ()4.131|) . 
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t=-4 



t=-5 




Figure 9.7: Snapshots of S;^ and v in a diagonal G2 cosmology, showing 
step-like structures in v. 
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Figure 9.8: Snapshot of the variables (E;^, E^, v) in a diagonal G2 cosmology 
at t = -4. 
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9.3.2 Spiky structures 

We now illustrate spiky structures in the trigger 
and provide evidence to support the claim that the 
by (mH): 

A = fC a.ic{T^Qi). 
Recall that OT G2 cosmologies satisfy ()5.56p : 

S2 = 0. (9.10) 

Consider the following initial data at t = with 7 = | (in conjunction 
with (jnH) and ^^): 

e = 0.1, (^1^)0 = 1, ^ = 0, (S^)o = 0.2, (Sx)o = 0.2, 

(iV_)o = 0, ^]o = 0.6, (fiA)o = 0.1, S2 = 0. (9.11) 

We ran CLAWPACK with 512 grid points, from t = to t = —10, storing 
data at intervals of t = 0.1. Figure . 1 II shows the snapshots of S_ and the 
trigger at t = —2, t = —6 and t = —8. We see the strong suggestion that 
spikes develop as t —00. Figures 15 . 1 21 shows the snapshots of the variables 
{T,^,T,^ , N^) in the (T,^,!^^, N^) space. Each dot represents the state at 
a grid point. We see that the solution shadows the Taub vacuum Bianchi II 
orbits on the sphere ()6.70p : 

1 = (iV^)2 + (S^)2 + (S^)2 

along some grid points, as it approaches the stable Kasner arc between T3 
and Qi. This may be compared with the Taub vacuum Bianchi II orbits 
shadowed by the transformed WM solution in Figure 16.21 There are two 
points stuck on the arc to the left of the T3 Taub Kasner point, creating two 
spikes. Along all other grid points, the orbits approach the past attractor 
(19. 9|) . The points closest to T3 take the longest to traverse the Taub vacuum 
Bianchi II orbits. 

The simulations also confirm that 

lim aiX = 0, (9.12) 

t— >— 00 



in OT G2 cosmologies, 
past attractor is given 

(9.9) 



along spike points, despite large spatial gradients. Thus asymptotic silence 
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t=-2 



t=-6 



t=-8 



0<x<2ji 0<x<2ji 



Figure 9.11: Snapshots of S_ and the trigger at t = —2, t = —6 and 
t = —8, showing the asymptotic signature (S_)sig(a;) and the formation of 
permanent spikes in A^_. 



holds along spike points. Our result has been mentioned in Andersson et al. 
2004 page S50]. 
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Figure 9.12: Snapshots of the variables (E^, E?, N^) in the (E^, E?, N^) 
space. Each dot represents the solution at a grid point. 
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9.4 Effects of spatial inhomogeneity II: role 
of spatial derivative terms 

In this section, we illustrate the role of spatial derivative terms in creating 
spike transitions and shock waves. 



9.4.1 Spike transitions 

Consider the following initial data at t 

my- 



= with 7 = I (in conjunction with 



e = 0.1, (Si^)o = l, Ao = 0, (S_)o = 0.2, (Sx)o = 0.2, 

(iV_)o = 0, fio = 0.6, (fiA)o = 0.1, (S2)o = 0.1. (9.13) 

High resolution is essential to simulate a spike transition correctly. We 
ran CLAWPACK with 16384 grid points, from t = to t = —4, storing data at 
intervals of t = 0.1. ^ 

Figure 19.131 shows the orbits along three grid points: two typical points 
at grid number 3000 and 3500, and the spike point at grid number 3267. The 
typical points undergo a sequence of three transitions, while the spike point 
undergoes a spike transition. The spike transition orbit is a straight line, 
taking the solution to the Kasner point that is reached by the sequence of 
A^_-Ex-A^- transitions, in about the same time. 

Figure 19.141 shows the snapshots of H^, the triggers , and the 
spatial derivative terms (Ei^)^ d^N^ , (Ei^)^ dx'^^ on the interval between 
grid point 3000 and grid point 3500. A spike forms at t = —2.6 in all of 
the variables, but also activates at the spike point (grid point 3267), and 
develops a very steep gradient. In addition, (Ei^)^ dxN^ and (Ei^)^ dx^^ 
develop spikes that have bigger amplitudes than the spikes in both and 
S^. A reversal of the spike in occurs some time between t = —2.6 and 
t = —2.7. The dynamics appears to be inherently inhomogeneous. 

We now take a closer look at the dynamics. Figures l9.15H9.l'ni show 
the /^-normalized terms in the evolution equations ()D.6|) - ()D.7|) for the pair 
(Ex,A^-). First, a spike forms in with the same mechanism as perma- 
nent spikes. From Figure I9.15t we see that diN^ is significant enough 



^The simulation in Section I^TI does not have high enough resolution to correctly sim- 
ulate spike transitions, but was acceptable for the typical background SH dynamics. 



CHAPTER 9. NUMERICAL SIMULATIONS 



152 




Figure 9.13: The orbit of an N_ spike transition along Xgpike; from t = —2.2 
to t = —3.3, projected onto the (S^, T,^) plane, shown as a solid line. The 
orbits of a sequence of A^_-Sx-A^- transitions along two nearby grid points 
are shown dashed line and a dash- dot line. 

that, when added to the algebraic terms on the right hand side of the 9jSx 
equation, it leads to a sign change in dtT^x at the spike point, which then ac- 
tivates Sx at the spike point to start the spike transition. From Figure 1^.1 61 
we see that an active Ex at the spike point in turn gives a ^iSx large enough 
that, when added to the algebraic terms on the right hand side of the dtN_ 
equation, it leads to a sign change in dtN^ at the spike point, which causes 
the reversal of the spike in A^_. 

What separates transient spikes like these from the permanent spikes that 
we observed in OT G2 cosmologies? Here, diN^ is large enough to activate 
Sx at the spike point. In permanent spikes, diN_ is too small, so Sx still 
decays to zero at the spike point. Whether diN_ is large enough depends on 
the Kasner point where the spike occurs (see Section [731) • 

The terminology spike transition is introduced by the author and was 
first referred to by Andersson et al. 2004 [31 page S50], who gave a survey of 
regular transitions in vacuum OT G2 cosmologies, using Hubble-normalized 
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Figure 9.14: An iV_ spike transition. Hubble-normalized variables are plot- 
ted. A reversal of the spike occurs between t — —2.6 and t — —2.7. 
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Figure 9.15: Influence of Ei^d^N^ on dtJ^x- At t = -2.4, Ei^d^N^ is large 
enough along the spike point that dtT,^ changes sign, causing Ex to develop 
a spike. 
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Figure 9.16: Influence of Ei^d^T^x on dtN_. The Ex spike that was started 
by a large Ei^dxN_ t — —2.4, now becomes large enough aX t — —2.6 
that the slope of dtN_ changes sign, causing the reversal of the N_ spike at 
t = -2.7. 
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variables. By examining the interplay between and Ex , we shed new light 
on the mechanism of spike transitions, building on the work of Garfinkle & 
Weaver 2003 |25j . What we have seen here is the transition of a so-called high- 
velocity spike. ^ In OT G2 cosmologies, a high-velocity spike is by definition 
a transient spike, which either reappears as a permanent, low- velocity spike 
eventually or disappears 

Similarly, a transient spike can occur in Ex, and will undergo a spike 
transition that is equivalent to a sequence of Sx-A^--Ex transitions. 

The simulation also suggests that spike transitions may violate asymptotic 
silence: as t — > — 00, tends to zero, but r becomes large during a spike 
transition, due to the term A^-Sx (see equation ()D.14|l ). ^ In other words, 
during a spike transition, the orbit does not shadow the silent boundary. If 
spike transitions occur indefinitely as t — —00, then asymptotic silence is 
violated. Inadequate numerical spatial resolution confines the simulations 
to very short times, so we have not been able to accurately simulate two 
consecutive spike transitions. ^ We shall leave this for future research. 

There are still questions about what happens to low-velocity spikes in 
generic G2 cosmologies. Will they be permanent or will they undergo tran- 
sitions? We have been unable to specify initial conditions that produce such 
transitions before numerical resolution becomes inadequate. We shall also 
leave this for future research. 



^The so-called velocity is equal to VSy S?. + S-^j ^'^d a spike transition is called a 

"slide-down" in the velocity in because + S-^ decreases by a predicted amount 
during a spike transition. Garfinkle & Weaver also recognized that a sign change in the 
evolution equations is crucial for a spike transition. 

^When doing linearization on the Kasner circles, we dropped the spatial derivative 
terms. The linearized evolution equations for and Sx indicate that N_'E,y^ tends to 
zero into the past. This linearization does not hold when the spatial derivative terms are 
large enough, as in the case of spike transitions. 

'^For a fixed grid size Ax, the largest gradient of a jump in a variable Y that can be 
represented is (Ymax — ^min) / Ax. A finer grid size is needed if the numerical gradient of a 
spike is observed to approach this upper bound and i?i^(ymax — ymin)/Ax is observed to 
be growing, otherwise the numerical solution can underestimate diY and thus can delay 
or even prevent a spike transition from occurring. This can mislead one to predict that 
spike transitions will eventually stop as t ^ —00, or that diY tends to zero as i — > — 00 
(see Andersson et al. 2004b Figure 2). 
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9.4.2 Shock waves 

Numerical experiments with shock waves show that the constraints (Cc)3 and 
(Cfs) (see ()D.16|) - ()D.17|) ) are unstable numerically near shock waves. To avoid 
this problem, we shall set fix and S2 to zero. As mentioned in Section IH21 
the presence of diA on the right hand side may prevent the system from being 
symmetric hyperbolic, which may be bad for simulations of shock waves. To 
avoid this problem we shall set A = 0. 

Consider the following initial data at t = with 7 = | (using the scheme 
for the case A = in Section 18.11 

A = = S2 = 
Ei^ = 100, S_=0.1, Sx=0.01, 
A^_ = 0.01, A^x=0.01, fi = 0.01, = 0.5 + 0.3sina; . (9.14) 

Here, a large Ei^ is chosen to create inhomogeneities with a short- wavelength 
(relative to the horizon). The amplitude in the fluctuation of v is also chosen 
to be large. These conditions are favourable for the formation of shock waves. 

We ran CLAWPACK with 512 grid points, from t = to t = —1, and from 
t = to t = 1 storing data at intervals of t = 0.01. Figures 17H9. 18l show 
the formation of a shock wave in the tilt variable v into the past and into 
the future respectively. 

The shock wave subsequently smoothes out or dissipates asymptotically. 
The mechanism of this dissipation is not clear. SH dynamics cannot dissipate 
a discontinuity, but would a small ^iX term dissipate it? 

Figure ID. 191 shows a snapshot of all the variables into the future (at t = 
0.06). We see that shock waves also appear in the profiles of fl, E+, q and r, 
but not in Sx, N— and Nx-, which seem to become non-differentiable at the 
location of the shock wave, but still continuous. The profiles of Ei^ and S_ 
seem smooth, but it could be that their non-differentiability is barely visible. 
Also plotted in Figure 19.191 are the /^-normalized Weyl scalar S'^ — Ti^ and 
the Hubble-normalized Weyl scalar {S^Y — {H^Y- 



CHAPTER 9. NUMERICAL SIMULATIONS 158 




0.5 


1 

0.5 


1 

0.5 


1 

0.5 




t=-0.02 



t=-0.1 



t=-0.5 



t=-0.9 



0.5 


1 

0.5 


1 

0.5 


1 

0.5 



t=-0.04 



t=-0.2 



t=-0.6 



< X < 27C 



< X < 2jt 



< X < 271 



t=-1 



< X < 27C 



Figure 9.17: Shock wave in v into the past. 
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9.18: Shock wave in v into the future. 
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Y=4/3, t=0.06 
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Figure 9.19: Snapshot of all variables, showing shock waves or non- 
differentiability into the future. 
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9.5 Future asymptotics 

In this section, we illustrate the future asymptotic behaviour of the tilt vari- 
able V, as predicted by ()5.7()|) . 

Consider three sets of the following initial data at t = 0, with 7 = 1.2, | 
and 1.5 (in conjunction with (j9.1|) ): 

e = 0.09, (^1^)0 = 1, Ao = 0.1, (S_)o = 0.2, (Sx)o = 0.2, 
(iV_)o = 0.2, fio = 0.4, (^]A)o = 0.01, (£2)0 = 0.1. (9.15) 

Changing 7 affects only the initial condition for v. 

We ran ICN with 512 grid points, from t = to t = 10, storing data at 
intervals of t = 0.1. 

As predicted by ()5.70p . the limit of the tilt variable v depends on the 
parameter 7. Figure 19.201 plots the snapshots of v for the three solutions. 
This confirms that if 7 < |, f tends to zero; if 7 = |, f tends to a function of 
X that depends on the initial condition; for 7 > |, f tends to ±1, depending 
on the initial condition. * 

As mentioned in Section 19.11 the evolution is asymptotically silent as 
t 00: 

Ei^ ^0, r , (9.16) 

and the solutions tend to 

nA = i, (^],A,s+,s_,Sx,S2,A^-,iVx) = 0. (9.17) 

The numerical decay rates agree with the analytical results of Lim et al. 
2004 j^. Table gives the predicted growth rates from |1H], and the typical 
absolute errors of the numerical decay rates. The errors are small, indicating 
excellent agreement between numerical and analytical results. Note that 
for the case |, since v tends to an arbitrary function, we take f (10, x) as the 
signature function, and determine the decay rate of f (t, x)—v{lQ, x) dXt = 5. 



^That step-like structures can develop in the case 7 > | was first predicted by Rendall 
2003 jffll. 
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Figure 9.20: Snapshots of v for three solutions with 7 
showing the asymptotic signature Vsig{x) in the case 7 > 
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Table 9.2: Predicted and actual growth rates ai t — 10. Errors are obtained 
along typical timelines, ignoring irregular growth rates along timelines where 
a variable crosses zero. 



Variable 


predicted growth rate 


Absolute error 


(order) 






7 = 1.2 




7 = 1.5 


E^' 


-1 


10-^ 


10-^ 


io-« 


A 


-1 


10-^ 


10-^ 


10-^ 


r 


-1 


10-5 


10-3 


10-3 


s+ 


-2 


10 


10 


10-3 


A/_ 


—2 


10-5 


10-5 


10 




-2 


10-3 


10-3 


10-3 


S2 


-3 


10-9 


10-9 


10-9 


N_ 


-1 


10-5 


10-5 


10-5 




-1 


10-5 


10-5 


10-5 




-2 


10-3 


10-5 


10-3 


n 


-37 
-4 
-4 


10-5 


10-5 


10-3 


V 


37-4 


10-3 






V - Vsig{x) 


-1 

2(37-4) 
2-7 




10-2 


10-^ 



9.6 Close-to-FL epoch 

In this section, we describe the evolution of a cosmological model that is 
close to the flat FL solution at t = 0. By "close to flat FL" we mean the 
Hubble-normalized shear and the Weyl scalars are sufficiently small, and Q 
is sufficiently close to 1 (see WE, page 62). 

Consider the following initial data a,t t — for an OT G2 cosmology 
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(S2 = 0) with 7 = 1 (using the scheme for the case A = in Section IHTTj) . 

Ei^ = l, S_ = O.OlsinSx, Sx = 0.01 sin 5x, = 0, A^x = 0, 

V = 0.01 sinx, n = 0.98G+, = 1 + (7 - l)v'^, 

r = -2(0.98)(0.01)sina;, 
= 10"^exp[-;^(0.98)(0.01) sin^ |]. (9.18) 

The initial value for Q is chosen so that the initial value for is close to 
zero (see Figure IU.24j) . 

We ran CLAWPACK with 512 grid points, from t = to t = —15, storing 
data at intervals of t = 0.1, and from t = to t = 15, storing data at intervals 
of t = 0.005. 

Figure 19.211 shows the evolution along the grid point 100 of 512. The 
graphs of In(ii^i^)^ and show three distinct epochs - two asymptotic and 
one intermediate. They are the Kasner asymptotic epoch from t = —15 to 
t = —4, the flat FL intermediate epoch from t = —4 to t = 3.48, and the de 
Sitter asymptotic epoch from t = 3.48 to t = 15. The fact that the fiat FL 
epoch is of finite duration is a reflection of the fact that the fiat FL model 
is unstable ^ into the future and into the past. The transition from the fiat 
FL epoch to the de Sitter asymptotic epoch is rapid, with steep slopes in 
the graphs of and Q^. The transition from the Kasner epoch to the flat 
FL epoch is less rapid, with a milder slope in the graph of . The growth 
rate of Ei^ is as expected, since Ei^ = (i?i^)oe^* for the Kasner solutions, 
Ei^ = (Ei^)oe* for the fiat FL solution with 7 = |, and = (^1^)06"* for 
the de Sitter solution. As t increases, Ei^ grows exponentially until ~ | 
at t = 3.48, when Ei^ reaches a maximum of 23.7902 and starts to decrease 
exponentially. 

Evolving into the past from t = 0, the solution enters the Kasner asymp- 
totic regime and approaches the Kasner arc {T3Q1). A snapshot of the spatial 
profile at t = —15 would show two spikes forming in Ex, and the asymptotic 
spatial signature, i.e. the dependence of S_ on x. This asymptotic behaviour 



^In the Hubble-normalized SH state space, the flat FL solution is described by an 
equihbrium point that is a saddle point, reflecting its instability within the SH context. 
The simulations also reflect unstable inhomogeneous modes. 

-"^"Notice that the {Cp)^ constraint has an unacceptably large error at about t = 3.48. As 
a result, the simulation loses accuracy temporarily. As mentioned in Section [8.4l stability 
of constraints is still an open issue in numerical simulations. 




Figure 9.21: Evolution along a fixed timeline - from the past asymptote to 
the future asymptote. 
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is similar to the one shown in Figure IH. Ill 

Evolving far into the future, the solution approaches the de Sitter solution. 
A snapshot of the spatial profile at t = 15 would show the asymptotic spatial 
signature, i.e. the dependence of v on x. This asymptotic behaviour is 
qualitatively the same as the one shown in Figure I9.2(JI 

We now take a detailed look at the close-to-FL intermediate epoch. Evolv- 
ing into the future from t = 0, we observe oscillations in the variables. Fig- 
ure 121221 shows in detail the evolution during the close-to-FL epoch and the 
earlier part of the de Sitter epoch, during which these oscillations occur. 
We see 45 oscillations in the pair (Sx,A^-), 27 in A^x), and 5 in v. 
How do these numbers relate to the wavelengths of spatial inhomogeneities? 
Figure 19.251 shows the snapshot at t = 3.48, revealing that the wavelength 
of spatial inhomogeneities is 27r/5 in the pair (Ex,iV_), 27r/3 in (E_,A^x), 
and 2n in {Q,v). (S_, iVx, Sx, iV-) are gravitational field variables, whose 
characteristic speed of propagation is the speed of light, while f is a mat- 
ter field variable, whose characteristic speed of propagation is the speed of 
sound, given by 1/7 — 1 times the speed of light. Comparing the numbers 
of oscillations with the wavelengths and taking into account the speeds of 
propagation, we have excellent agreement - we see | as many oscillations in 
A^x) as in (Sx, A^-), because the wavelength of spatial inhomogeneities 
in A^x) is | times as long. Since 7 = |, we can predict that there should 

be ^1 — 1 X I 0.5774 x | as many oscillations in v as in (E_,A^x); be- 
cause waves in v travel 0.5774 as fast as A^x) do, and the wavelength of 
spatial inhomogeneities in f is 3 times as long. The simulations agree with 
the prediction. 

We also observe that the oscillations are most rapid aXt = 3.48, when Ei^ 
reaches its maximum. This is expected, because the frequency is equal to 
the speed of propagation divided by the wavelength. Recall from ()5.85|) that 
the coordinate speed of light is Ei^Af. Since we have chosen a t-coordinate 
such that = 1, for a fixed wavelength the frequency is proportional to £"1^. 
We note that the frequency depends on the choice of t-coordinate, and the 
physical frequency should be computed with respect to the clock time, for 



-'^-'^ There are also oscillations in O, approximately of amplitude 0.02, but these are not 
visible in the graph because of the scale. 

^^For characteristic eigenfields and their corresponding characteristic velocities of prop- 
agation, see van Elst et al. 2002 77. 




Figure 9.22: Evolution along a fixed timeline - details of the close-to-FL 
epoch. 
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which M = p. Since (5 ^ oo into the past and /? ^ y |- into the future, the 
physical frequency actually decreases with time. 

With the shear and spatial curvature variables staying close to zero during 
the flat FL epoch and subsequently during the de Sitter epoch, the instability 
of the fiat FL model appears to be primarily due to the growth of Ei^ into the 
future, leading to growth in the Weyl curvature, which describes anisotropy in 
the free gravitational field. The spatial derivative terms are the largest at 
t = 3.48, when Ei^ attains a maximum. The largest spatial derivative terms 
are Ei^d^^x and Ei^dxN_, both of which have magnitudes approximately 
equal to 



magnitude ~ max ^i?i^(Sx , A^-) j x 



27r 

wavelength 
(24 X 0.005) X 5 = 0.6 . 



Thus the Weyl scalars £ and H (see ()4.117p and ()5.35|) - ()5.44|) ). which contain 
spatial derivative terms, are approximately | x 0.6 = 0.2 at t = 3.48. This 
growth in the Weyl curvature scalars mirrors the growth of Ei^ and the 
decrease of the wavelength of inhomogeneities relative to the particle horizon. 
From t = 3.48 onwards, the Weyl curvature scalars decrease, reflecting the 
decrease of Ei ^ and the increase of the wavelength of inhomogeneities relative 
to the event horizon. Figure 19.231 shows the evolution of the Weyl scalars 
along the grid point 100 of 512, confirming the above argument. 

We have performed a variety of simulations of G2 cosmologies with a 
close-to-FL epoch. Decreasing the initial size of Q\ will increase the dura- 
tion of the fiat FL epoch, defined hj Q ^ 1, S ^ 0, into the future. On the 
other hand, this change will allow more time for Ei^ to grow, which will in 
turn increase the Weyl scalars significantly. Increasing the initial value of Ei^ 
will also increase the Weyl scalars. We comment that there are three modes 
(with two fields in each mode) that develop oscillations, namely {Q, v) de- 
scribing acoustic waves, and (S_, A^x) and (Sx, A^-) describing gravitational 
waves. In the initial condition ()9.18|) we introduced three different wave- 
lengths into the three modes (sinx, sin 3a; and sin 5a;). Coupling between 
these three modes does occur, but is rather weak. For example if f = ini- 



^■^Such behaviour can also occur in SH cosmologies of Bianchi type VIIq (see Nilsson et 
al. 1999 'ST'), i.e. the Hubble-normalized Weyl curvature can be large even though the 
Hubble-normalized shear is close to zero. 
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Figure 9.23: Evolution of the Hubble-normalized Weyl scalars along a fixed 
timeline, showing growth as Ei^ increases. 

tially the gravitational wave modes will subsequently turn v on, with sin 3x 
and sin 5x modes. We note that the two fields in each mode can have different 
wavelengths. 

In conclusion we note that the simulations during the close-to-FL epoch 
—4 < t < 3.48 can be viewed as non-linear perturbations of the flat FL model. 
It would be interesting to use simulations of this nature to test the validity 
of the linear perturbations of the flat FL model that are used extensively in 
cosmology. 
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Figure 9.24: Snapshot of the initial condition ai t — 0, showing a small E+. 
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7^4/3, t=3.48 




0<x<2jc 0<x<2jc 



Figure 9.25: Snapshot at t = 3.48, showing the wavelengths in the different 
modes. 



Chapter 10 
Conclusion 



In this thesis, by constructing a unified framework, analyzing explicit solu- 
tions, providing heuristic arguments, and performing numerical simulations, 
we have made considerable progress in understanding the dynamics of inho- 
mogeneous cosmologies with one spatial degree of freedom, the so-called G2 
cosmologies. We have highlighted the role of spatially homogeneous cosmolo- 
gies, in particular the tilted models, in providing the background dynamics. 

The dynamics of tilted SH cosmologies is still not well-understood. In 
Chapter El we have given a unified classification scheme for G2 cosmologies 
and for G'2-compatible SH cosmologies, based on the G2 action and on the 
tilt degrees of freedom. In Chapter 0J this classification leads to a unified 
system of equations and spatial gauge choice that is valid for G'2-compatible 
SH cosmologies of all Bianchi types. By analyzing the stability of the Kasner 
equilibrium points, we have shown that the occurrence of Mixmaster dynam- 
ics is determined by the G2 action, and that the tilt determines whether the 
fiuid vorticity is dynamically significant near the initial singularity. 

The notion of asymptotic silence provides the link between G2 and SH 
dynamics. This link enabled us to establish that Mixmaster dynamics occurs 
in a G2 cosmology if and only if the G2 action is generic, by making use 
of the analysis of Mixmaster dynamics in SH cosmologies. We have also 
given the past attractors for subclasses of G2 cosmologies that do not permit 
Mixmaster dynamics. 

We have made significant progress in understanding the development of 
steep spatial gradients. In Chapter El we have given a detailed analysis of 
two classes of explicit solutions, the first one with step-like structures, and 
the second one with spikes. The second one is new, and is the only known 
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explicit G2 solution with spikes. We have made a detailed comparison of the 
asymptotic size of different terms in the evolution equations, which revealed 
that spatial derivative terms are not always negligible. Furthermore, the 
growth of spatial derivative terms triggers the so-called spike transitions. As 
a result, the BKL conjecture II does not hold in general for G2 cosmologies. 

All numerical simulations of G2 cosmologies to date have considered vac- 
uum models, and have used a metric formulation. Chapter El is one of two 
pieces of work, the other being Andersson et al. 2004 (4j [2j , in which numeri- 
cal simulations of G2 cosmologies have been performed in a dynamical system 
framework. Our work complements that of Andersson et al. in that they con- 
sider the past asymptotic dynamics of vacuum G2 cosmologies, whereas we 
have considered the full range of the dynamics for non-vacuum models in 
which the matter content is a radiation perfect fluid and a positive cosmo- 
logical constant. 

Our simulations show that, for models that are close to FL at some time, 
the evolution consists of three distinct epochs, namely the past and future 
asymptotic epochs and the close-to-FL intermediate epoch. The numerical 
results support the following conclusions, within these epochs. 

• Past asymptotic silence, future asymptotic silence (if A > 0), the BKL 
conjecture I, and the Kasner Attractivity Conjecture are valid along 
almost all orbits of a typical G2 cosmology. 

• Mixmaster dynamics occurs along a typical timeline of a typical G2 cos- 
mology in the past asymptotic regime, provided that the group action 
is generic. 

• The dynamics in asymptotically silent regimes consists of the back- 
ground SH dynamics and a spectrum of inhomogeneous phenomena, 
in which the trigger variables and the wavelengths of inhomogeneities 
play crucial roles. 

• A trigger variable that has both signs causes step-like or spiky struc- 
tures to form. 

• A spike transition occurs if the spatial derivative of a trigger causes the 
time derivative in the evolution equation to change sign. 

• Sufficiently short wavelength spatial inhomogeneities in the matter vari- 
ables can lead to the formation of shock waves. 
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• A long enough closc-to-FL epoch allows the wavelength of spatial in- 
homogeneities to decrease, which leads to wave-like phenomena and 
potentially a large Weyl curvature. 

• Wave-like phenomena disappear as the solution enters the past asymp- 
totic regime or the de Sitter future asymptotic regime, which is asymp- 
totically silent. 

We conclude by commenting on future research that this thesis leads to. 
Firstly, the unified formulation that we have given provides a framework 
for proving the conjectures, in particular the asymptotic silence conjecture, 
the BKL conjecture I, and the Kasner Attractivity Conjecture, that appear 
to govern the asymptotic dynamics of tilted SH and G2 cosmologies. Sec- 
ondly, the analysis of the tilted SH cosmologies can now be completed using 
the unified formulation. Thirdly, it is possible that further analysis of the 
close-to-FL simulations, which represent non-linear perturbations of the flat 
FL model, may shed light on an important question: how reliable are the 
standard linear perturbation analyses of the flat FL model, which play a 
fundamental role in cosmology? 



Appendix A 
Go action 



In this appendix, we prove Proposition 13.11 in Section 13.11 about G2 models 
with a hypersurface-orthogonal KVF. We spht the proposition into two parts: 
Proposition lA.il and its corollary. 

Proposition A.l. Consider a G2 model presented in a group-invariant orbit- 
aligned frame with 62 and 63 tangent to the G2 orbits. Then the G2 admits 
a hypersurface-orthogonal KVF if and only if 



a^^n^^ = 0, nABaiV^ = 0, ^a^^a^V^ = 0. (A.l) 
Proof : Let ^ be a HO KVF, i.e. 




0. 



(A.2) 



Expanding this yields 



l^oi^c = 0. 



(A.3) 
(A.4) 
(A.5) 



Using ()2.42|) . equations ()A.5|) simplifies to 



aiA^^ = 0. 



(A.6) 
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The group-invariant frame yields 

eo(e'')+7Ve^ = (A.7) 
ei(e^)+7Ve^ = 0, (A.8) 

which are then contracted with ^c^"^ b and substituted into ()A.3|1 and ()A.4jl 
to obtain 

*aABeV = and nABi^i^ = Q (A.9) 
respectively. Equations ()A.6|) and ()A.9|) combine to yield ()A.1|) . □ 

"<^=": Assuming ()A.1|) . use the rotation freedom to set ^ 

n_ = 0. (A.12) 

Then dHJ yields 

CX12 = 0, a^= 0. (A.13) 

The evolution equations for ?7,_, au and Ux, and the (6*0)2 constraint then 
imply that 

n+ = 0, ^]l = 0. (A.14) 

Equations ()A.12|) - ()A.14|) imply that the frame vector 62 is hypersurface- 
orthogonal. We also have, trivially, 

f^i = —VSax, n+ = VSn^. (A. 15) 

Using Lemma 1 below, this means 62 is parallel to a KVF. Thus this KVF is 
also hypersurface-orthogonal, and is independent of the choice ()A.12|) . □ 



^For convenience we introduce the following notation: 



2n+, hAB = V3("- , (A.IO) 



and 

(712 = \/3cr3, (713 = \/3ct2. (A.U) 
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Lemma 1 

Consider a G2 model presented in a group-invariant orbit- aligned frame with 
62 and 63 tangent to the G2 orbits. Then 

fii = —VSax, n+ = -\/3n_ (A. 16) 

if and only if 62 is parallel to a KVF. 

Proof. Let ^ and r} be KVFs. As seen in the proof of the Proposition, the 
group-invariant frame yields 

eo(e'')+7Ve^ = (A.17) 

ei(e'')+7Ve^ = 0, (A.18) 

and hkewise for rj. 

"<^": Assuming ^ is parallel to ea, we have ^3 = 0. Then (|XT7|l - (|08|l 
imply that 

7^02 = = 7^2, (A.19) 

or equivalently ()A.16|1 . as required to prove. □ 

Assuming ()A.16|) . or equivalently, 

7^02 = = 7^12, (A.20) 

we obtain the following: 

eote) + 7'o3e3 = 0, 61(^3) + 7^13^3 = 0, (A.21) 

and likewise for rf. Equation ()A.21|) has a very special form. If ^3 = at 
one point in spacetime, then ()A.21|) implies that ^3 = at every point in 
spacetime. 

If ^3 = or ?73 = 0, then ^ or 7/ is parallel to e2 and we are done. If not, 
pick any point P and a constant K so that 

e3 + i^r/3 = at P. (A.22) 

But ^3 + Krj^ satisfies the equations 

eo(e3+i^r73)+7'o3fe+^^3) = 0, ei(e3+i^r73)+7'i3(e3+i^r73) = 0, (A.23) 



APPENDIX A. G2 ACTION 



178 



which now imphes 

^3 + K7]s = (A.24) 
at every point in spacetime. Thus ^+-^'57 is parallel to 62 and we are done. □ 

Corollary of Proposition lA.ll 

Consider a G2 model presented in a group-invariant orbit-aligned frame with 
62 and es tangent to the G2 orbits. If the G2 admits a hypersurf ace- orthogonal 
KVF, then 

e^V^iB = 0, e^%At^iB = 0, *^^Bt^^V' = 0. (A.25) 



Proof. Using ()A.6|1 and ()A.9|) . we simplify ()A.7|) and ()A.8|) to the following: 

eo(e^) = n^eAB^'' (A.26) 



ei(e^) = WeAB^- (A.27) 



We shall then differentiate (TOl) and (TOl) . First, using and (TOHIl . 

eolc^iA^^) = yields 

TTi^e^ = 0, (A.28) 

which implies 

e^'^T^iAaiB = 0. (A.29) 
Next, using (jTK^ and (fOTjl . ei(aiA^^) = yields 

g^e^ = 0, (A.30) 

which implies 

e^^'qACTiB = 0. (A.31) 

Similarly, using (E3n|l and (jS^ll), eo(*aAB^^^^) = and ei(n^BC^^^) = 
combine to yield 

"'f^ABi'-e = 0, (A.32) 

which implies 

*^ABa^^(y^^ = 0, (A.33) 
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as required. It also follows from (jA.lj) and ()A.25j) that 

*frABa^'' = 0, nABn^"" = 0. ^ (A.34) 

□ 



^To show this, use the notation IjA.lOp and introduce the notation 

(n_,nx) = (n cos (j),n sin (f)), (1712, ctis) = (cr cos 6*, a sin 6*), 
and similarly for ttab- 



Appendix B 
Rotation formulae 



For G2 cosmologies, the spatial gauge freedom is represented by a position- 
dependent rotation of the spatial frame vectors {62,63}: 

:::)(::). ^-^(*-)- ^ 

For G'2-compatible SH cosmologies, the rotation is restricted to = 

Writing the commutators using ()B.1|) gives the formulae for the transfor- 
mation of the variables (regardless of /?- or Hubble-normalization) under a 
rotation: 

R 

(I) 

it) 
it) 

Other variables, namely Ei^, A, E+, r, q, Q, vi and ^a, are not affected by 
the rotation. 



N+ + di<f) 



cost/) 

— sin 

COS0 

— sin 

cos 20 

— sin 20 



sm</) 

COS0 

sin0 

COS0 

sin ; 
cos! 



cos 20 sin 20 
— sin 20 cos 20 



S3 

V3 



(B.2) 
(B.3) 

(B.4) 
(B.5) 
(B.6) 
(B.7) 
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Note that the determinant of Nab, given by 

det Nab = NI - 3(7V! + K), (B.8) 

is invariant if = (f){t), but will change if (p also depends on x. The sign of 
det Nab is thus significant for G2-compatible SH cosmologies, but not for G2 
cosmologies. 



Appendix C 

Asymptotic expansions at spike 
points 



For a general vacuum OT G2 solution, Kiclienassamy & Rendall 1998 
equations (5)-(6)] gave a generic expansion for the metric functions P and 
Q (see 

P = -2k{x)t - (t)[x) + ■ ■ ■ (C.l) 
Q = Xo(x) + e^^'^^Hi^x) + ■ • • , (C.2) 

where k{x) satisfies < k{x) < 1. The above expansion is inadequate 
in satisfying the evolution equations ()6.34|) - ()6.37|) . To obtain an adequate 
expansion, we start with the vacuum Bianchi II solution 



where F = Foe"^^^"*^-**. We replace the constants k and Fq by arbitrary 
functions of x and integrate the evolution equations ()6.34|) - ()6.37|) to obtain 
updated N^, Sx, -/V_ and E_ in succession. Finally, we integrate the asymp- 
totic expansion for E_ and Ex to obtain the asymptotic expansion for P and 

Q, using dnaoi). 
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The necessary expansion for P and Q is as follows: 

P = -2k{x)t - 0(x) + A(x)e''=(")* + B{x)e^^^-^^' + C(x)e«'=(")* + ■ ■ ■ (C.4) 
Q = Xo{x) + ^(x)e^'=(")* + ^(x)e«'=(")* + x(x)e^2^(")* 

+ e(x)te^* + C(a;)e'* + ■ ■ ■ , (C.5) 



where 



e = -Ai,, x = A'i,, s, = -^^y (C.7) 
^ - + 16(1^ " 8(1^ " 4(1^^ • ^^-^^ 



Rendall & Weaver 2001 [021 transform the Kichenassamy-Rendall expan- 
sion ()C.1|) - ()C.2|) to produce expansions for P, Qt and for the transformed 
solution.^ We shall transform ()(I4j) - ()(I5|l to obtain a consistent expansion 
for the transformed solutions (S_, A^x , Sx , A^-)- The resulting expansions 
for the transformed solution are as follows. Along x = a;spike(— oo), where 
Aro(a;spikc(-oo)) = 0: ^ 

E_ = -^(1 + k) + ^(1 - A;)5e^(i-^)* + ^(e^^fct ^ ^sii-k)t _ ^^u^ ^^^g^ 
iVx = -^A;' (-te^*) + -^(2X^^e-2^ + ^Oe^* + 0{-te^^^^^''^') (C.IO) 
= ^.i.X^e-V^^-'^)* + 0(e2(i+'^)* + e^^i-^)*) (C.ll) 
iV_ = -|(l-A;)ee-^(-te2(2-'=)*) 

- ^e-<^ [e + 4(1 - A;)C + 8(1 - k)Bi^] e^^^-'^)* 



^The temporal and spatial coordinates used by Rendall & Weaver 2001 |^ are r and 
xrwi which are related to ours according to r = — 2i, xjiw — 2x. For example, Xq(x) 
here would correspond to 2Xq{xrw) in jf>2| . 

^It turns out that the asymptotic expansions for 9tS_, dtN^, dt^y, and dtN_ are the 
same as the time derivative of the asymptotic expansions (jC.9p - (|C.12|) . and so we do not 
list them. 
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(C.13) 



__i=A;'e2* + - k) [X'^ie'^'t' + 2k' B] {-te^i^-^"^') 

- Ts^o^e-^^ + ^(1 - k)B' - ^[4(1 - k)X',^|Je-'^ + k']B 



_ ^(1 _ k)X'^Ce^^f 62(3-2^=)* 



(C.14) 
(C.15) 



(C.16) 



Along typical points: ^ 



X. 



0^0 



(C.17) 
(C.18) 



iV_ 



1 YV-'/'^2(l-fc)t 



2(l+fc)t _^ g6{l-fc)t 



+ C(-te 



(C.19) 
(C.20) 



The above expansion is consistent with the evolution equations ()6.34|) - ()6.37p . 
To simplify to the expansion for the transformed WM solution in Section 17^ 
set 

k = \, = 0, V' = /i, ^o = -2/i, ^0=4/2. (C.21) 

In this case, the A{x) and B{x) terms have the same power and they cancel 
each other in the expansion for P, since A{x) = —B{x) from ()C.6|) and 

The above expansion is a new result, and provides the asymptotic expan- 
sion in the framework of the orthonormal frame formalism. 



■^It turns out that the asymptotic expansions for the temporal and spatial derivatives of 
E_ , Nx , S X and are the same as the temporal and spatial derivatives of the asymptotic 
expansions (|(y.l7p - (|(].20p . and so we do not list them. 



Appendix D 



Equations for numerical 
simulations of G2 models 

Following the derivation in Section 15.11 we arrive at a system of equations 
()5.8p - ()5.29p with /^-normalized variables for G2 cosmologies. As mentioned 
at the beginning of Chapter |S1 we shall restrict our considerations to G2 
cosmologies with one tilt degree of freedom. In doing so, the shear spatial 
gauge ()5.46p simplifies to ()8.2p : 

S3 = 0, = V3A^_, R = -VsSx . (D.l) 

We drop the index on the remaining tilt variable Vi, and simplify the evolu- 
tion equation ()5.33|) for v'^ by collecting the di{v^), di InG^ and diVi terms 
into an Ei^dxV term. As explained in Chapter |H1 we solve the constraint 
(I5.2ip to determine S4.. 

Since Ei^ and A have the same evolution equation, we can save numerical 
resources by writing 

A = B{x)Ei^ (D.2) 

and evolve Ei^ only, and obtain B{x) from the initial conditions. See Sec- 
tioning for more discussion of the initial conditions. 

Numerical experiments show that S2, and v are prone to getting 
out of range, e.g. f2 becoming negative. To prevent this we evolve lni?i\ 
lnE2, lnf2 and arctanh v instead. 

The resulting system for numerical simulations, given below, consists of 9 
evolution equations ()D.3|) - ()D.11|) and two constraint equation ()D.16|) - ()D.17|) . 
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Evolution equations: 

dJnEi^ = g + 3E_ 



(D.3) 



dtE_ + Ei^d^Ny, = {q + 3E+ - 2)E_ + 2ANy, + 273(^2 - iv!) - VsSs^ 

(D.4) 

dtN^ + ^i^a.S. = {q + 3S+)iVx (D.5) 
^fSx - Ei^d^N_ = {q + 3S+ - 2 - 2y3S_)Sx - 2AA^_ - 2^A^xiV- 

(D.6) 

dtN^ - Ei^d^Ey, = {q + 3S+ + 2V3E_)iV_ + 2V3ExiVx (D.7) 
^tlnSs = g- 2 + y3S_ (D.8) 



7G_(l-t;2)^, 



G. 
27 



£^1 9j;(arctanh v) 



+ 



+ 1) - i(l - 3S+)(1 + _ 1 + yl^; 



(D.9) 



9t(arctanh v) + ^ — ^ — ^Ei^d^^lnQ 



-(7G_)-i [(2 - 7)G+r - 7^(87 - 4 + 3(2 - 7)E+] 



+ 27(7 - 1)^^;'] 



dtlnQA = 2{q + 1) , 



(D.IO) 
(D.ll) 
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where 



2(1 - A^) + Ei'd.A - 3E2' - 1(2 - 7) 



-fl-SflA (D.12) 



l-A'^ + lEi^d^A + A^n + 2A{N^j:_-N_J:^) 



G. 



n 



Gi = 1 ± (7 - ly 

Constraint equations: 



/[2{1-A')] (D.13) 
(D.14) 
(D.15) 



= (Cch = ^/<9,S2 -{r + 3A- VSN^)^^ 
= (Cf,) = {E^d, - 2r)0A . 



(D.16) 
(D.17) 



For simulations in the dc Sitter asymptotic regime (where I^a — 1 and q-\-l 
tend to zero), we use VL^ — 1 and g + 1 as variables to maintain numerical 
precision. 
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